Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Prior for R_0 -- suggestion #105

Open
maxbiostat opened this issue May 13, 2020 · 5 comments
Open

Prior for R_0 -- suggestion #105

maxbiostat opened this issue May 13, 2020 · 5 comments

Comments

@maxbiostat
Copy link

Hi, first of all, many many thanks for the awesome work all of you have been doing!

Secondly, I have a minor question:

What's the justification for the prior currently employed for R_0? The hierarchy

kappa ~ half-normal(0, 0.5)
R_0 ~ half-normal(m, kappa)

leads to a "witch hat"-looking prior on R_0.

image

Would a simpler Gamma/log-normal prior not be better here? Easier to elicit and probably captures key features such as Pr(R_0 > 1) more easily. If you're interested I can help construct such a prior and then submit a PR, but I have no easy way of running the model to check results as of now.
@flaxter
Copy link
Collaborator

flaxter commented May 13, 2020 via email

@flaxter
Copy link
Collaborator

flaxter commented May 13, 2020 via email

@flaxter
Copy link
Collaborator

flaxter commented May 13, 2020 via email

@maxbiostat
Copy link
Author

maxbiostat commented May 13, 2020 •
edited
Loading

Hi @flaxter ,

It's more likely I'm the buggy one:

Code

LogMean <- function(realMean, realSD){
  ## takes REAL mean and REAL standard deviation; returns LOG mean
  realVar <- realSD^2
  mu <- log(realMean/ sqrt(1 + (realVar/realMean^2)))
  return(mu)
}
LogVar <- function(realMean, realSD){
  ## takes REAL mean and REAL standard deviation; returns LOG variance
  realVar <- realSD^2
  sigmasq <- log(1 + (realVar/realMean^2))
  return(sigmasq)
}
###################
M <- 1e6
kappa_ic <- abs(rnorm(n = M, mean = 0, sd = 1/2))
R0_ic <- abs(rnorm(n = M, mean = 3.28, sd = kappa_ic))

kappa_bram <- abs(rnorm(n = M, mean = 0, sd = 1))
R0_bram <- abs(rnorm(n = M, mean = 2.4, sd = kappa_bram))

m <- 3 # E[R0]
v <- .9^2 # Var(R0)

# Gamma
alpha <- m^2/v
beta <- m/v

#Log-normal
mu <- LogMean(realMean = m, realSD = sqrt(v) )
sigma <- sqrt( LogVar(realMean = m, realSD = sqrt(v) )    ) 

R0_gamma <- rgamma(n = M, shape = alpha, rate = beta)
R0_ln <- rlnorm(n = M, meanlog = mu, sdlog = sigma)

forplot <- data.frame(R_0 = c(R0_ic,
                              R0_bram,
                              R0_gamma,
                              R0_ln),
                      study = rep(c("Imperial",
                                    "BRAM-COD",
                                    "Gamma",
                                    "Log-normal"), each = M))

library(ggplot2)

ggplot(forplot, aes(x = R_0, colour = study, fill = study)) +
  geom_density(alpha = .4) +
  scale_x_continuous(expression(R_0), expand = c(0, 0)) +
  scale_y_continuous("Density", expand = c(0, 0)) +
  geom_vline(xintercept = 1, linetype = "longdash") +
  theme_bw(base_size = 16)

summy <- function(x, na.rm = TRUE){
  ans <- data.frame(mean = mean(x, na.rm = na.rm),
                    sd = sd(x, na.rm = na.rm),
                    lwr = quantile(x, probs = .025, na.rm = na.rm),
                    uprr = quantile(x, probs = .975, na.rm = na.rm))
  return(ans)
} 


lapply(list(bramcod = R0_bram, imperial = R0_ic, Gamma = R0_gamma, LN = R0_ln), summy)

Last line gives

$bramcod
         mean       sd       lwr      upr
2.5% 2.435088 0.911255 0.5271803 4.580456

$imperial
         mean        sd      lwr      upr
2.5% 3.280113 0.4989465 2.190974 4.374173

$Gamma
         mean        sd      lwr      upr
2.5% 2.999272 0.8991958 1.503212 4.995466

$LN
        mean        sd      lwr      upr
2.5% 3.00032 0.8997761 1.617084 5.106452

BRAM-COD is a recent Brazilian study with a model based on your model from report 13.

@maxbiostat
Copy link
Author

I realised the script above has a (as of yet inconsequential) bug in the generation of the random samples. Technically, the model specified in Stan is a truncated normal, not a folded normal, so the correct generation of the samples would be something like

M <- 1e6
library(truncnorm)
kappa_ic <- rtruncnorm(n = M, mean = 0, sd = 1/2, a = 0, b = Inf)
R0_ic <- rtruncnorm(n = M, mean = 3.28, sd = kappa_ic, a = 0, b = Inf)

kappa_bram <- rtruncnorm(n = M, mean = 0, sd = 1, a = 0, b = Inf)
R0_bram <- rtruncnorm(n = M, mean = 2.4, sd = kappa_bram, a = 0, b = Inf)

Giving

> lapply(list(bramcod = R0_bram, imperial = R0_ic, Gamma = R0_gamma, LN = R0_ln), summy)
$bramcod
         mean        sd       lwr    uprr
2.5% 2.471883 0.9082468 0.6756234 4.67399

$imperial
         mean       sd      lwr     uprr
2.5% 3.280978 0.497174 2.192721 4.371233

$Gamma
         mean        sd      lwr     uprr
2.5% 2.999731 0.9005028 1.502073 5.006691

$LN
        mean        sd      lwr     uprr
2.5% 3.00048 0.9005886 1.616915 5.110658

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
2 participants
@flaxter @maxbiostat