Bayes rule for data analysis:

\[\underbrace{P(\theta \, | \, D)}_{posterior} \propto \underbrace{P(\theta)}_{prior} \times \underbrace{P(D \, | \, \theta)}_{likelihood}\]

Markov Chain Monte Carlo

• sequence of representative samples from \(X \sim P\)

JAGS

• declarative language to specify data-generating model
  • model = prior + likelihood function
  • samples of whatever variable we want:
    • conditioning on data \(\rightarrow\) additional constraints on log-score

• fun with JAGS models
  • hierarchical models / latent mixture models
  • think: "data-generating process"
  • ideas & notation
• priors in hierarchical models
• shrinkage
• examples, examples, examples

require('rjags')
modelString = "
model{
  theta ~ dbeta(1,1)
  k ~ dbinom(theta, N)
}"
# prepare for JAGS
dataList = list(k = 7, N = 24)
# set up and run model
jagsModel = jags.model(file = textConnection(modelString), 
                       data = dataList,
                       n.chains = 2)
update(jagsModel, n.iter = 5000)
codaSamples = coda.samples(jagsModel, 
                           variable.names = c("theta"),
                           n.iter = 5000)

ms = ggs(codaSamples)
ggs_density(ms) + 
  stat_function(fun = function(x) dbeta(x, 8,18), color = "black") 

• rational Bayesian agent
• originally believes \(\theta \sim Beta(1,1)\) (any bias equally likely)
• observes \(k = 7\) successes in \(N = 24\) tosses

model{
  theta ~ dbeta(1,1)
  k ~ dbinom(theta, N)
}

• rational Bayesian agent
• originally believes \(\theta \sim Beta(16,16)\)
  • e.g., has observed 15 successes and 15 failures
  • has updated believes via Bayes rule
• observes another \(k = 7\) successes in \(N = 24\) tosses

model{
  theta ~ dbeta(6,16)
  k ~ dbinom(theta, N)
}

• boundedly-rational, Bayesian agent with imprecise memory
• has observed 15 successes and 15 failures but doesn't clearly remember
• observes another \(k = 7\) successes in \(N = 24\) tosses

model{
  a ~ pois(5)
  b ~ pois(15)
  theta ~ dbeta(a+1,b+1)
  k ~ dbinom(theta, N)
}

• latent variables can depend on other latent variables
• likelihood can be a (direct) function of only some latent variables

dummy

\[ \begin{align*} P(\theta, a, b \, | \, D) & \propto P(\theta, a, b) \ P(D \, | \, \theta, a, b)\\ & = P(a) \ P(b) \ P(\theta \mid a, b) \ P(D \, | \, \theta) \end{align*} \]

dummy dummy

\[ \begin{align*} P(\theta, s \, | \, D) & \propto P(\theta, s) \ P(D \, | \, \theta, s)\\ & = P(s) \ \underbrace{P(\theta \mid s)}_{\text{degenerate!}} \ P(D \, | \, \theta) \end{align*} \]

model{
  s ~ dbern(0.5)
  theta = ifelse(s == 1, 0.75, 0.25)
  k ~ dbinom(theta, N)
}

• type of variable:
  • continuous: circle
  • discrete: square
• information status:
  • observed: grey shade
  • latent: white
• dependency:
  • stochastic: single line
  • deterministic: double line

• priors are complex objects
• should make sense at every level in the hierarchy:
  • check dependent priors of intermediate latent variables!
• can have interesting effects
  • e.g., shrinkage (later in this class)

model{
  b ~ dgamma(1,1)
  a ~ dgamma(1,1)
  theta ~ dbeta(a,b)
  k ~ dbinom(theta, N)
}

dummy dummy

  ggplot(data.frame(x = c(0.0001,10)), aes(x)) +
         stat_function(fun = function(x) dgamma(x,1,1))

• \(\text{Beta}(a,b)\) has mode \(\omega\) and concentration \(\kappa\): \[ \begin{align*} \omega & = \frac{a-1}{a+b-2} & \kappa &= a + b \end{align*} \]
• reformulate: \[ \begin{align*} a & = \omega(\kappa - 2) + 1 & b & = (1-\omega)(\kappa-2)+1 \end{align*} \]

model{
  omega ~ dbeta(1,1)
  kappaMin2 ~ dgamma(0.01,0.01)
  theta ~ dbeta(omega*kappaMin2+1, 
                (1-omega)*kappaMin2+1)
  k ~ dbinom(theta, N)
  kappa = kappaMin2+2
}

dummy dummy

• hierarchical models with latent parameters:
  • individual-level differences
  • item-level differences
  • …

dummy dummy

• TT practitioners sense energy fields
• they sense which hand is near another persons body

TTdata = readr::read_csv('../data/02_therapeutic_touchers.csv') %>% 
         group_by(s) %>% 
         summarize(k = sum(y), N = length(y))
# k - number of successes per subject
# N - number of trials per subject
head(TTdata)

## # A tibble: 6 × 3
##       s     k     N
##   <chr> <int> <int>
## 1   S01     1    10
## 2   S02     2    10
## 3   S03     3    10
## 4   S04     3    10
## 5   S05     3    10
## 6   S06     3    10

ggplot(TTdata, aes(x  = s, y = k/N)) + geom_bar(stat='identity')

ggs_density(ms, family = "omega|kappa")

msTheta = filter(ms, ! Parameter %in% c("omega", "kappa"))
msTheta$subject = with(msTheta,
                       factor(regmatches(Parameter, regexpr("[0-9]+", Parameter)), 
                              levels = 1:28))  
ggplot(msTheta, aes(x = subject, y = value)) + geom_boxplot() + ylab("bias level (posterior)")

• careful: "shrinkage" is a misleading term!
  • shrinkage of variance in low-level parameters (compared to non-hierarchical)
  • variance can also increase, depending on the particular model
• hierarchical modeling "binds together" low-level parameters
  • this is an additional assumption (can be good or bad)
• hierarchical "binding" reduces "freedom" of free variables
  • can guard us against overfitting
• hierarchical group-level binding gives us population-level estimates
  • generalize to new observations

Candidates from Thailand and Moldovia take part in a quiz. Questions are about the history of each country. Here's the data of who got which question right.

show(CountryQuizData)

##    QA QB QC QD QE QF QG QH
## P1  1  0  0  1  1  0  0  1
## P2  1  0  0  1  1  0  0  1
## P3  0  1  1  0  0  1  0  0
## P4  0  1  1  0  0  1  1  0
## P5  1  0  0  1  1  0  0  1
## P6  0  0  0  1  1  0  0  1
## P7  0  1  0  0  0  1  1  0
## P8  0  1  1  1  0  1  1  0

model{
  alpha ~ dunif(0,1)    # chance correct if origin matches question
  beta ~ dunif(0,alpha) # mismatch   
  for (i in 1:dim(k)[1]){
    x[i] ~ dbern(0.5)
  }
  for (j in 1:dim(k)[2]){
    z[j] ~ dbern(0.5)
  }   
  for (i in 1:nx){
    for (j in 1:nz){
      theta[i,j] = ifelse(x[i] == z[i], alpha, beta)
      k[i,j] ~ dbern(theta[i,j])
    }
  }   
}

• hierarchical models: chains of dependencies btw. variables
• information from data percolates "upward":
  • stronger impact on nearby nodes
  • less impact on higher nodes (convergence checks!)
• high-level nodes can unify estimates for lower-level variables
  • sometimes called "shrinkage"
  • better conceived as: group-level prior homogeneity assumption
• complex data generating models:
  • prior + likelihood
  • our assumptions about how the data could have been generated

Friday

• boot camp on complex models in JAGS

Tuesday

• model comparison

obligatory

• skim Lee & Wagenmakers chapters 5 & 6
  • decide what you want to explor in class
• finish homework 3