Fundamentals of Bayesian Data Analysis in R

Course bayesian methods definition: A method of figuring out unobservable quantities given known facts that uses probability to describe the uncertainty over what the values of the unknown quantities could be.

*4 hours, 23 Videos, 58 Exercises

Strengths of bayesian methods:

• Learns from data
• Used for hyp. testing and linear regression.
• Flexible enough to handle problem-specific models.
  • can include information sources in addition to the data.
  • can make comparisons between groups or data sets.
  • can use the result of a bayesian analysis to do decision analysis
  • can change the underlying statistical model.

Prior = a probability distribution that represents what the model knowns before seeing the data.

Posterior = a probability distribution that represents what the model known after having seen the data.

Bayesian Inference requires data, priors, and a generative model.

• a generative model is a mathematical expression with a set of rules used to feed parameter values and simulate data.

Example: Probability that we have found a vaccine that cures the zombie outbreak.

• using the binomial process as a generative model

# parameters
p_success <- .15 # probability of success
n_zombies <- 13 # trails

# simulate data
data <- c()
for(zombie in 1:n_zombies){
  # save if prop is between 0 & 1 and greater than our prob. of success
  data[zombie] <- runif(1, min=0, max=1) < p_success 
}

data <- as.numeric(data)
data
mean(data) # our posterior is close to our theorized probability of success of .15
sum(data) # number of successes

binomial function

• Parameters - probability and trial size.

N=13

# prob. of success, close to our prior belief
sum(rbinom(n=1, size=N, prob=.15)) / N 

Lots of simulations

set.seed(909)

# prob. of success, close to our prior belief
sum(rbinom(n=100000, size=100, prob=.10)) / (100000*100) 
hist(rbinom(n=100000, size=100, prob=.10))

Adding an uninformed prior distribution, includes out uncertainty in a binomial distribution.

• Ex. we want to know how many times our add will be clicked on, clicks=success, we are not sure if it is 10% of the time so we will assign a prior that is equally likely as low as 0% and as high as 20%

# number of observations
n_samples <- 100000 

# number of trails
n_ads_shown <- 100 

 # probability of success with a uniform prior representing our uncertainty
proportion_clicks <- runif(n = n_samples, min = 0.0, max = 0.2)

n_visitors <- rbinom(n = n_samples, size = n_ads_shown, prob = proportion_clicks)

# posterior probability, close to our previous prob constrained at .10
sum(n_visitors) / (100000*100) 

# probability of visitors, not a uniform dist like when  probability was constrained to .10
hist(n_visitors) 

# a uniform distribution of our uncertainty
hist(proportion_clicks) 

# probability of 5 clicks or more
sum(as.numeric(n_visitors >= 5)) / n_samples 

Updating our prior:

Bayesian inference = conditioning on data, in order to learn about parameter values.

• What would be the unknown parameter value conditional on 7 visitors? .47
• What would be the unknown parameter value conditional on 11 visitors? .87

prior <- data.frame(n_samples, n_ads_shown, proportion_clicks, n_visitors)
head(prior)

Now we collected data on 100 observations and 13 clicked our add, we can update our prior to calculate the probability of getting 5 clicks or more.

posterior  <- prior[prior$n_visitors == 13, ]

# updated posterior
head(posterior) 

prior <- posterior

head(prior)

n_samples <-  nrow(prior)

n_ads_shown <- 100

prior$n_visitors <- rbinom(n=n_samples, size = n_ads_shown, 
                           prob=prior$proportion_clicks)

hist(prior$n_visitors)

# new probability of 5 clicks or more after updating our belief from our posterior
sum(prior$n_visitors >= 5) / length(prior$n_visitors)  

Note. This approach would scale bad if we had more data or complexity.

Using the beta distribution

• Takes 2 parameters, alpha and beta.

beta_sample <- rbeta(n = 1000000, shape1 = 1, shape2 = 1)
hist(beta_sample)
head(beta_sample)

beta_sample <- rbeta(n = 1000000, shape1 = 10, shape2 = 10)
hist(beta_sample)

# changing the shape of the beta distribution

in our add example if most adds get clicked 5% of the time, but some adds get clicked as low as 2% and as high as 8%, what beta parameters would best capture this prior distribution?

hist(rbeta(n=100, shape1 = 5, shape2 = 95))

New tweak to our previous model, adding an informed prior and data to inform the new posterior

n_draws <- 100000
n_ads_shown <- 100

# Change the prior on proportion_clicks
proportion_clicks <- rbeta(n_draws, shape1 = 5, shape2 = 95) # NEW CHANGE
n_visitors <- rbinom(n_draws, size = n_ads_shown, prob = proportion_clicks)
prior <- data.frame(proportion_clicks, n_visitors)
posterior <- prior[prior$n_visitors == 13, ]

# This plots the prior and the posterior in the same plot
par(mfcol = c(2, 1))
hist(prior$proportion_clicks, xlim = c(0, 0.25))
hist(posterior$proportion_clicks, xlim = c(0, 0.25))

Comparing between groups or datasets; example video (13 out of 100 clicked) vs text (6 out of 100 clicked) ads

n_draws <- 100000
n_ads_shown <- 100
proportion_clicks <- runif(n_draws, min = 0.0, max = 0.2)
n_visitors <- rbinom(n = n_draws, size = n_ads_shown, 
                     prob = proportion_clicks)
prior <- data.frame(proportion_clicks, n_visitors)

# Create the posteriors for video and text ads
posterior_video <- prior[prior$n_visitors == 13, ] # 13 of 100 clicked
posterior_text <- prior[prior$n_visitors == 6, ] # 6 of 100 clicked


# Visualize the posteriors
hist(posterior_video$proportion_clicks, xlim = c(0, 0.25))

hist(posterior_text$proportion_clicks, xlim = c(0, 0.25))

Compare difference in posterior proportions between video and text ads

posterior <- data.frame(
    video_prop = posterior_video$proportion_clicks[1:4000],
    text_prop  = posterior_text$proportion_click[1:4000])
    
# Calculate the posterior difference: video_prop - text_prop
posterior$prop_diff <- posterior$video_prop - posterior$text_prop

# Visualize prop_diff
hist(posterior$prop_diff)

# Summarize prop_diff
summary(posterior$prop_diff)

Most likely difference, and probability that the proportion of clicks us larger for the video than text ad.

median(posterior$prop_diff)
sum(posterior$prop_diff > 0)/length(posterior$prop_diff) 

Note. the posterior mean or median is the summary distribution over some parameter, while the mean or median of a data is a summary of the data.

Credible interval, contains the underlying parameter with a certain probability.

Decision analysis: taking a result of a statistical analysis and post process to make an informed decision (i.e., save lives, make more money).

• example: video ads cost 25 cents, text 5 cents, and we make on average $2.53 per visitor.

video_cost <- .25
text_cost <- .05
visitor_spend <- 2.53

# P(people  that click the video ad * how much we made on visitors - cost)
posterior$video_profit <- posterior$video_prop * visitor_spend - video_cost 

# P(people that click the text ad * how much we made on visitors - cost)
posterior$text_profit <- posterior$text_prop * visitor_spend - text_cost 

# difference in profits this is the distribution we care about
posterior$profit_diff <- posterior$video_profit - posterior$text_profit 
head(posterior)
hist(posterior$profit_diff)

# Calculate a "best guess" for the difference in profits
median(posterior$profit_diff)

# Calculate the probability that text ads are better than video ads
sum(posterior$profit_diff < 0) / length(posterior$profit_diff)

In looking at the histogram there is no real difference in profits over lost yet text ads are marginally cheaper to run.

Changing the underlying statistical model

• binomial will not work here.
• instead, we use the poisson distribution modeling the number of clicks per day.
• Explore the poisson distribution. Example: on a cloudy day what is the probability of breaking even if we sell on average 11.5 ice creams?

ice_cream <- rpois(n=10000, lambda = 11.5) # lambda is our average count
hist(ice_cream)
sum(ice_cream >= 15)/ length(ice_cream)

unlikely that we will break even.

Applying it to our ad problem.

• We now get charged by the day, rather than by ad, to run an ad. Our clicks are now per day instead of per ad.
• new day shows that in 1 day we received 19 clicks.
• how many daily clicks should we expect on average?

n_draws <- 100000
n_ads_shown <- 100

# flat prior from 0 to 80 clicks per day
mean_clicks <- runif(n_draws, min = 0, max = 80) 

# 2 parameters size and mean
n_visitors <- rpois(n_draws, lambda = mean_clicks) 
                     
# probability of getting 19 clicks
posterior <- prior[prior$n_visitors == 19, ] 

summary(posterior$mean_clicks)

quantile(posterior$mean_clicks, prob=c(.05, .95))

Probability rules

• Any number between 0 and 1
• A statement (un)certainty ##### Probability of drawing one of the four aces in a stack of 52 cards

prob_to_draw_aces <- 4/52
prob_to_draw_aces

Probability of drawing all four aces in a row

prob_to_draw_four_aces <- (4/52) * (3/51) * (2/50) * (1/49)
prob_to_draw_four_aces

Calculating probabilities with functions rather than simulating data, conditioning on 10%

• conditional probability of getting 13 clicks on our adds given that the success rate is 10%

n_visitors <- rbinom(n = 100000, size = 100, prob = 0.1)
sum(n_visitors == 13) / length(n_visitors)

# or 
dbinom(13, 100, .10)

probability of getting more than 13 clicks

n_visitors <- rbinom(n = 100000, size = 100, prob = 0.1)
sum(n_visitors > 13) / length(n_visitors)

# or 
pbinom(13, 100, .10, lower.tail=FALSE)

Conditioning on the data

n_ads_shown <- 100
n_visitors <- seq(0, 100, by=1)
proportion_clicks <- seq(0, 1, by=.01)

pars <- expand.grid(proportion_clicks=proportion_clicks, 
                    n_visitors=n_visitors)

pars$prior <- dunif(pars$proportion_clicks, min=0, max=.2)

pars$likelihood <- dbinom(pars$n_visitors, size=n_ads_shown, 
                          prob=pars$proportion_clicks)

pars$probability <- pars$likelihood * pars$prior

pars$probability <- pars$probability / sum(pars$probability)
sum(pars$probability)
head(pars)

# conditioning on data when visitors equal 6
pars <- pars[pars$n_visitors == 6,]

# normalize again
pars$probability <- pars$probability / sum(pars$probability)


# Plot the posterior pars$probability
plot(pars$proportion_clicks, y=pars$probability, type="h")

Bayes theorem:

\[ \underbrace{p(\theta|D)}_{posterior} ~=~ \frac{\overbrace{p(D|\theta)}^{likelihood}~ \overbrace{p(\theta)}^{prior}} {\underbrace{p(D)}_{evidence}} \]

The probability of different parameters \(\theta\) given some data \(D\) \(=\) The likelihood: the relative probability of the data given different parameter values \(P(D|\theta))\) \(\times\) The prior: the probability of different parameters before seeing the data \(P(\theta)\) \(\div\) The total sum of the likelihood weighted by the prior \(\sum{P(D|\theta)\times P(\theta)}\)

Mathematical notation for our model, convenient way to define model.

\[ \eta_{ads} = 100 \]

\[ P_{clicks} \sim Uniform(0.0,0.2) \]

\[ \eta_{visitors} \sim Binomial(\eta_{ads}, p_{clicks}) \]

for the poisson model

\[ \mu_{clicks} \sim Uniform(min:0, max:80) \]

\[ n_{visitors} \sim Poisson(\mu_{clicks}) \]

working with the normal distribution

• 2 parameters \(\mu\) and \(\sigma\)

# mean of 20, sd of 2
rnorm(n=5, mean=20, sd=2) 

# temperature readings from a lake
temp <- c(19, 23, 20, 17, 23) 

# relative likelihood of our data points (not probabilities)
likelihood <- dnorm(x= temp, mean=20 , sd=2) 
likelihood

# multiply all numbers together (small number)
prod(likelihood) 

# working on the log scale to avoid poor scability of data
log(likelihood)

example: how smart are zombies? Data are from IQs of zombies.

• What we’re interested in is how much we can learn about the mean zombie IQ from this data.

\[ \mu \sim Normal(mean:100, sd:100) \]

\[ \sigma \sim Uniform(min:0, max:50) \]

\[ IQ_i \sim Normal(\mu, \sigma) \]

\[ iq = 55, 44, 34, 18, 51, 40, 40, 49, 48, 46 \]

fit the model

# The IQ of a bunch of zombies
iq <- c(55, 44, 34, 18, 51, 40, 40, 49, 48, 46)

# Defining the parameter grid
pars <- expand.grid(mu = seq(0, 150, length.out = 100), 
                    sigma = seq(0.1, 50, length.out = 100))

# Defining and calculating the prior density for each parameter combination
pars$mu_prior <- dnorm(pars$mu, mean = 100, sd = 100)
pars$sigma_prior <- dunif(pars$sigma, min = 0.1, max = 50)
pars$prior <- pars$mu_prior * pars$sigma_prior

# Calculating the likelihood for each parameter combination
for(i in 1:nrow(pars)) {
  likelihoods <- dnorm(iq, pars$mu[i], pars$sigma[i])
  pars$likelihood[i] <- prod(likelihoods)
}
pars$probability <- pars$likelihood * pars$prior
pars$probability <- pars$probability / sum(pars$probability)

levelplot(probability ~ mu*sigma, data=pars)

mean is roughly around 42 with a sd of + or -20

The posterior distribution for our example

• pars contains the data frame representing the posterior zombie IQ.

head(pars)
sample_indices <- sample( 1:nrow(pars), size = 10000,
    replace = TRUE, prob = pars$probability)
head(sample_indices)

# Sample from pars to calculate some new measures
pars_sample <- pars[sample_indices, c("mu", "sigma")]

# Visualize pars_sample
hist(pars_sample$mu)

# Calculate the 0.025, 0.5 and 0.975 quantiles of pars_sample$mu

quantile(pars_sample$mu, c(0.025, 0.5, .975))

We estimate the mean zombie IQ to be 42 (95% CI: [35, 50])

what range of zombie IQs should we expect? And how likely is it that the next zombie you encounter is, at least, moderately intelligent?

pred_iq <- rnorm(10000, mean = pars_sample$mu, 
                 sd = pars_sample$sigma)
# Visualize pred_iq
hist(pred_iq)

# Calculate the probability of a zombie being "smart" (+60 IQ)
pred_iq <- rnorm(10000, mean=pars_sample$mu, sd=pars_sample$sigma)

# the Pr that the next zombie you'll meet will have an IQ of >=60 
sum(pred_iq >= 60)/length(pred_iq) 

very unlikely that the next zombie we encounter will have an IQ of 60 or more.

BEST: package

• assumption, data comes from a t-distribution and uses MCMC
• 3 parameters, mean, mu, and degree of freedom (small df assumes more outliers)

fit <- BESTmcmc(iq)
fit

We get a close mean and sd from previous analysis.

• the following plot looks similar to our previous histogram

plot(fit)

Comparing 2 data sets.

• example: zombies that have a regular diet and those that have a brains based diet

iq_brains <- c(44, 52, 42, 66, 53, 42, 55, 57, 56, 51)
iq_regular <- c(55, 44, 34, 18, 51, 40, 40, 49, 48, 46)
# Calculate the mean difference in IQ between the two groups
mean(iq_brains) -  mean(iq_regular)

# Fit the BEST model to the data from both groups
best_posterior <- BESTmcmc(iq_brains, iq_regular)
best_posterior

our mean difference was 9.3, our posterior mean difference (51.780 - 43.102)= 8.678

plot(best_posterior)

Bayesian estimate using a t-distribution is robust to outliers.

• suppose one of the regular diet zombies has a score of 175

iq_brains <- c(44, 52, 42, 66, 53, 42, 55, 57, 56, 51)
iq_regular <- c(55, 44, 34, 18, 51, 40, 40, 49, 48, 175)
# Calculate the mean difference in IQ between the two groups
mean(iq_brains) -  mean(iq_regular)

# Fit the BEST model to the data from both groups
best_posterior <- BESTmcmc(iq_brains, iq_regular)
best_posterior

t=There is more uncertainty in our estimate yet, we still conclude that brain eating zombies are more likely to have a higher IQ than regular diet zombies.

plot(best_posterior)