get_beta_params: Get parameters for Beta distribution

In gilesjohnr/propvacc: Functions to model and simulate vaccine dose coverage

Description

This function finds the two shape parameters for the Beta distribution of a random variable between 0 and 1. Note that the function uses a different method depending on the arguments supplied. The three methods are:

1. When the mean (mu) and variance (sigma) are supplied, the solution is found analytically.

2. When observed probabilities (probs) at each quantile (quantiles) are given, the solution is found by minimizing the Sum of the Squared Errors (SSE) using the Nelder-Mead optimization algorithm. Note that the fitting algorithm performs best when the five standard quantiles are supplied (Min, 25th, Median, 75th, Max).

3. When only observed probabilities (probs) are supplied, the function uses Maximum Likelihood Estimation (MLE).

Usage

get_beta_params(mu = NULL, sigma = NULL, quantiles = NULL, probs = NULL)

Arguments

mu	scalar giving the mean μ
sigma	scalar giving the variance σ^2
quantiles	vector of quantiles for which proportions are observed. Expects: c('min', '25th', 'median', '75th', 'max').
probs	vector of observed probabilities or proportions

Value

A list containing the two shape parameters of the Beta distribution

Author(s)

John Giles

See Also

Other beta_params: get_beta_params_age()

Examples

### Analytic solution
mu <- 0.6
sigma <- 0.005
prm <- get_beta_params(mu=mu, sigma=sigma)
curve(dbeta(x, prm$shape1, prm$shape2), 0, 1, lwd=2,
      xlab="Random variable", ylab='Probability density')
abline(v=mu, col='blue')
abline(v=c(mu-sqrt(sigma), mu+sqrt(sigma)), col='blue', lty=2)


### Numeric solution using SSE
quantiles <- c(0.001, 0.25, 0.5, 0.75, 0.999)
probs1 <- c(0.26, 0.57, 0.73, 0.81, 0.95) # WHO region AFR, <12mos, n=16
probs2 <- c(0.86, 0.86, 0.88, 0.92, 0.94) # WHO region AFR, >12mos, n=4

v <- rbind(probs1, probs2)
dimnames(v) <- list(Age=c('< 12mos', '> 12mos'),
                    Quantiles=c('min', '25th', 'median', '75th', 'max'))

par(mfrow=c(1,2))
prm <- get_beta_params(quantiles=quantiles, probs=probs1)
curve(dbeta(x, prm$shape1, prm$shape2), 0, 1, lwd=2,
      xlab="Vaccine effectiveness", ylab='Probability density')
abline(v=probs1, col='blue', lty=2)
mtext(c('min', '25th', 'median', '75th', 'max'), 3, at=probs1, line=0.5, col='blue', cex=0.75)
mtext('<12mos', 3, at=0.05, line=-1.5)

prm <- get_beta_params(quantiles=quantiles, probs=probs2)
curve(dbeta(x, prm$shape1, prm$shape2), 0, 1, lwd=2,
      xlab="Vaccine effectiveness", ylab='Probability density')
abline(v=probs2, col='blue', lty=2)
mtext('>12mos', 3, at=0.05, line=-1.5)


### Maximum Likelihood Estimation
get_beta_params(probs=rbeta(100, 1.5, 2))

# Measles vaccination coverage in Namibia
# WHO-UNICEF estiamtes from:
# http://apps.who.int/immunization_monitoring/globalsummary/estimates?c=NAM

mcv1 <- c(82, 80, 75, 85, 83, 82, 76, 74,
          75, 76, 73, 69, 63, 73, 70, 70,
          68, 58, 69, 65, 64, 59, 68)/100
mcv2 <- c(50, 32)/100

prm <- get_beta_params(probs=mcv1)
curve(dbeta(x, prm$shape1, prm$shape2), 0, 1)
abline(v=mcv1, lty=2, col='blue')

prm <- get_beta_params(probs=mcv2)
curve(dbeta(x, prm$shape1, prm$shape2), 0, 1)
abline(v=mcv2, lty=2, col='blue')

gilesjohnr/propvacc documentation built on Aug. 24, 2020, 3:20 a.m.