Please read this: Name your file homework3_lastname_firstname. This set is due 12/2 before class. The submission should be in Rmarkdown, and should go to Fabian with subject BDA: Homework 3 (submit both .Rmd and .pdf). We encourage discussing exercises with fellow students, but only individual solutions will be accepted.

1. Infering … yeah, what actually?

Check out the following code.

require('rjags')
require('reshape2')
require('ggmcmc')
set.seed(1638) 
coolData = data.frame(A = rnorm(100, mean = 100, sd = 10),
                      B = rnorm(100, mean = 110, sd = 15))
coolData = melt(coolData, variable.name = "Drug", value.name = "HeartRate")

# specify model
modelString = "
model{
  # priors
  muA ~ dunif(30,200)
  muB ~ dunif(30,200)
  varA ~ dunif(0, 1000)
  varB ~ dunif(0, 1000)
  tauA = 1/varA
  tauB = 1/varB

  # likelihood
  for (i in 1:200){
    HeartRate[i] ~ dnorm(ifelse(Drug[i] == 0, muA, muB), 
                         ifelse(Drug[i] == 0, tauA, tauB))
  }

  muDiff = muB - muA
  varDiff = varB - varA
}
"

# prepare for JAGS
dataList = list(Drug = ifelse(coolData$Drug == "A", 0, 1),
                HeartRate = coolData$HeartRate)

# set up and run model
jagsModel = jags.model(file = textConnection(modelString), data = dataList, n.chains = 3)
codaSamples = coda.samples(jagsModel, variable.names = c("muDiff", "varDiff"), n.iter = 10000)

  1. Describe what is going on in a few words, and try to guess the intention of this analysis. I.e., what is most likely the research question that this analysis is trying to answer?

  2. Run the analysis and report the result. Include (i) a summary of the sampling outcome, (ii) density plots of the relevant parameters (individually per chain), and (iii) checks on the quality of the sampling (trace plots, R-hat and effective sample size). Interpret your results (i.e., should we think that we have approximated the true posterior? is our sampling efficient?).

  3. Given your answer to the question above about the likely research question behind this analysis, what would be your conclusion given this analysis?

2. Simple linear regression

We will work with data about the speed of a car and the distance needed to stop.

# install.packages('car') # if you don't have this package yet
require('car')
## Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
## logical.return = TRUE, : there is no package called 'car'
head(cars) # this is the data set we work with
##   speed dist
## 1     4    2
## 2     4   10
## 3     7    4
## 4     7   22
## 5     8   16
## 6     9   10

It is intuitive that there should be some kind of functional dependency between the speed that a car travels at and the distance it takes for it to come to a halt when breaking. Indeed, we can verify this by plotting ‘distance’ against ‘speed’ and then adding a simple linear regression line, obtained here with the R function lm.

We could get this regression line directly from ggplot, but then we might not realize that this is actually a regression analysis. We also want to know the values of slope and intercept of the regression line later on.

require('ggplot2')
coeff = coef(lm(formula = dist ~ speed, data = cars)) # get slope & intercept of linear model
qplot(cars[,'speed'], cars[,'dist']) + xlab('speed') + ylab('distance') + theme_bw() + 
  geom_abline(intercept = coeff[1], slope = coeff[2], color = "skyblue")

We are going to write a simple Bayesian linear estimator for ‘distance’ as a function of ‘speed’. Concretely, we want to implement the following Bayesian model, where \(x_{i}\) is the \(i\)-th entry in the data frame for variable speed and \(y_{i}\) is the \(i\)-th entry for variable distance:

\[\beta_0 \sim \text{Norm}(0, 1000)\] \[\beta_1 \sim \text{Norm}(0, 1000)\] \[\sigma^2_{\epsilon} \sim \text{Unif}(0, 1000)\]

\[\mu_i = \beta_0 + \beta_1 x_i\] \[y_i \sim \text{Norm}(\mu_i, \sigma^2_{\epsilon})\]

To explain, \(\beta_1\) is the slope, \(\beta_0\) is the intercept. We assume that these give us a mean \(\mu_i\) for each \(x_i\). The value \(y_i\) is then assumed to be normally distributed around this \(\mu_i\).

  1. Implement this model in JAGS. Rember that JAGS wants precision not variance as the second argument to a normal distribution!

  2. Run the model with sufficient burn-in and samples. Report the usual diagnostics on the samples (summary, traceplots, R-hat, efficient sample size) to make sure that your analysis is sound.

  3. Plot the density of the posteriors of \(\beta_0\) and \(\beta_1\) and get the 95% HDIs for both variables. Compare your results with the coefficients returned by lm. (If you get completely different results, there’s something wrong.)

  4. Which information do you get from the Bayesian analysis that you do not get from calling lm?

3. Discrete therapeutic touchers

We introduced a latent-mixture model for the therapeutic touch data in class (following Kruschke). Here, we would like to explore a related model, which assumes that there are three types of touchers: (i) suckers, (ii) guessers, (iii) sensors. Suckers have a low chance of success, guessers operate at around chance level, and sensors have an above average chance. We would like to use the data to infer the latent touch-type of each subject.

More concretely, for each subject i we are going to infer a latent type type[i]. The prior for type type[i] is an unbiased categorical distribution (each type is equally likely): type[i] ~ dcat(c(1/3, 1/3, 1/3)). Touch-types come from three disinct populations, each with its own mode, for \(j \in \{1,2,3\}\) the type of toucher, omega[j] of a beta distribution for the success chance theta. Suckers have omega[1] = 0.25, guessers have omega[2] = 0.5 and sensors have omega[3] = 0.75. The type of a subject thus determines the mode of a beta distribution: if i is of type type[i] then the mode of the beta distribution from which i’s theta[i] is drawn is omega[type[i]]. We assume that there is a single concentration kappa that is shared among all three type-populations. Given mode (for each subject) and concentration parameter, we determine a theta[i] (for each subject) as in the previous touch-model from class. The theta[i] feeds into a binomial likelihood function as before.

  1. Load the from file TherapeuticTouchData.csv, which is part of DBDA2Eprograms.zip, which is in our Dropbox in folder Kruschke_2015. Reformat it just like we did in class.

  2. Draw the hierarchical model described above. To save time, feel free to draw by hand (legibly, please!), take a picture, and include it in your file.

  3. Implement the hierarchical model in JAGS. Run it (with sufficient burn-in to ensure convergence), and use the code below to plot your outcome:

msType = ms %>% filter(substr(x = Parameter, start = 1, stop = 4) == "type") %>%
    mutate(subject = factor(regmatches(Parameter, regexpr("[0-9]+", Parameter)), 
                              levels = 1:28))
ggplot(msType, aes(x = subject, fill = factor(value))) + geom_bar()
  1. What would you conclude from this analysis about whether therapeutic touch works? Do you think that this is a good analysis (why/why not)? (You could think about stuff like: shrinkage, over-fitting, model complexity, …)