# maximum-likelihood empirical bayes calculations
# for take-home test 2, problem 2(B)

# unbuffer the output

# here's a function that does the empirical bayes analysis

# y and V are the data vectors

# the function starts at sigma.hat = 0 and iterates
# until two successive values of sigma.hat differ
# by no more than a small positive number epsilon

# n is a vector containing the combined (treatment + control)
# sample sizes in the k = 6 studies

aspirin.case.study.empirical.bayes.calculations <- 
  function( y, V, epsilon, n ) {

  sigma.squared.hat.old <- 0

  W.hat.old <- 1 / ( V + sigma.squared.hat.old )

  mu.hat.old <- sum( W.hat.old * y ) / sum( W.hat.old )

  sigma.squared.hat.new <- sum( W.hat.old^2 * 
    ( ( y - mu.hat.old )^2 - V ) ) / sum( W.hat.old^2 )

  m <- 0

  while ( abs( sigma.squared.hat.new - 
    sigma.squared.hat.old ) > epsilon ) {

    sigma.squared.hat.old <- sigma.squared.hat.new

    W.hat.new <- 1 / ( V + sigma.squared.hat.new )

    mu.hat.new <- sum( W.hat.new * y ) / sum( W.hat.new )

    sigma.squared.hat.new <- sum( W.hat.new^2 * 
      ( ( y - mu.hat.new )^2 - V ) ) / sum( W.hat.new^2 )

    m <- m + 1

  }

  mu.hat <- mu.hat.new

  sigma.squared.hat <- sigma.squared.hat.new

  sigma.hat <- sqrt( sigma.squared.hat )

  W.hat <- W.hat.new

  W.hat.normalized <- W.hat / sum( W.hat )

  n.normalized <- n / sum( n )

  B.hat <- V / ( V + sigma.squared.hat )

  theta.hat <- ( 1 - B.hat ) * y + B.hat * mu.hat

  se.hat.mu.hat <- 1.0 / sqrt( sum( 1 / ( V + sigma.squared.hat ) ) )

  se.hat.theta.hat <- sqrt( V * ( 1 - B.hat ) )

  mle.results <- list( m = m, mu.hat = mu.hat, 
    se.hat.mu.hat = se.hat.mu.hat, sigma.squared.hat = 
    sigma.squared.hat, sigma.hat = sigma.hat, n = n, 
    n.normalized = n.normalized, W.hat = W.hat, 
    W.hat.normalized = W.hat.normalized,
    B.hat = B.hat, y = y, theta.hat = theta.hat, 
    se.hat.theta.hat = se.hat.theta.hat )
    
  return( mle.results )

}

# specify the data values

y <- c( 2.77, 2.50, 1.84, 2.56, 2.32, -1.15 )

V <- c( 1.65, 1.31, 2.34, 1.67, 1.98, 0.90 )^2

# pick a value for epsilon; 10^( -7 ) seems reasonable

epsilon <- 1e-7

# specify the treatment and control sample sizes n.T and n.C in the 
# k = 6 experiments, and calculate n = n.t + n.C

n.T <- c( 615, 758, 317, 832, 810, 2267 )

n.C <- c( 624, 771, 309, 850, 406, 2257 )

n <- n.C + n.T

# now call the function

print( empirical.bayes.results <- 
  aspirin.case.study.empirical.bayes.calculations( y, V, epsilon, n ) )

# m is the number of iterations needed to achieve convergence:

# $m

# [1] 35

# mu.hat is the MLE of mu (this should agree with what You got
# from 'optim')

# $mu.hat

# [1] 1.446869

# se.hat.mu.hat is the estimated standard error of mu.hat 
# (this differs from what You got from 'optim' in an interesting way)

# $se.hat.mu.hat

# [1] 0.8089829

# sigma.squared.hat is the MLE of sigma^2

# $sigma.squared.hat

# [1] 1.530753

# sigma.hat is the MLE of sigma (this should agree with what You got
# from 'optim')

# $sigma.hat

# [1] 1.237236

# n is a vector of combined sample sizes (treatment + control)
# in the k = 6 experiments

# $n

# [1] 1239 1529  626 1682 1216 4524

# n.normalized is n divided by the sum of the n vector; this shows
# what proportion of the overall number of patients in the entire
# meta-analysis was in each experiment

# $n.normalized

# [1] 0.11455251 0.14136464 0.05787722 0.15551036 0.11242604 0.41826923

# W.hat is defined in the test (see formula (20))

# $W.hat

# [1] 0.2351142 0.3079905 0.1427276 0.2315001 0.1834474 0.4272130

# W.hat.normalized is W.hat divided by the sum of the W.hat values;
# w.hat[ i ] represents the weight on observation y[ i ] in the
# weighted average of the y[ i ] values that defines mu.hat

# $W.hat.normalized

# [1] 0.15387125 0.20156544 0.09340856 0.15150600 0.12005779 0.27959095

# B.hat is the amount of weight given to y[ i ] in the weighted average
# defining theta.hat (see formula (23) in the test)

# $B.hat

# [1] 0.6400983 0.5285426 0.7815193 0.6456306 0.7191873 0.3460425

# y is the vector of mortality differences (control - treatment)
# for the k = 6 studies

# $y

# [1]  2.77  2.50  1.84  2.56  2.32 -1.15

# theta.hat is the MLE of the vector theta (see model (15) in the test)

# $theta.hat

# [1]  1.9230658  1.9433752  1.5327602  1.8413283  1.6920550 -0.2513731

# se.hat.theta.hat is an approximate standard error for theta.hat

# $se.hat.theta.hat

# [1] 0.9898648 0.8994821 1.0937609 0.9941332 1.0492369 0.7278087