Normal approximation for Bioassay model.

ggplot2, grid, and gridExtra are used for plotting, tidyr for manipulating data frames

library(ggplot2)
library(gridExtra)
library(tidyr)
library(MASS)

Bioassay data, (BDA3 page 86)

df1 <- data.frame(
  x = c(-0.86, -0.30, -0.05, 0.73),
  n = c(5, 5, 5, 5),
  y = c(0, 1, 3, 5)
)

Grid sampling for Bioassay model.

Compute the posterior density in a grid

  • usually should be computed in logarithms!
  • with alternative prior, check that range and spacing of A and B are sensible
A = seq(-1.5, 7, length.out = 100)
B = seq(-5, 35, length.out = 100)
# make vectors that contain all pairwise combinations of A and B
cA <- rep(A, each = length(B))
cB <- rep(B, length(A))
# a helper function to calculate the log likelihood
logl <- function(df, a, b)
  df['y']*(a + b*df['x']) - df['n']*log1p(exp(a + b*df['x']))
# calculate likelihoods: apply logl function for each observation
# ie. each row of data frame of x, n and y
p <- apply(df1, 1, logl, cA, cB) %>% rowSums() %>% exp()

Sample from the grid (with replacement)

nsamp <- 1000
samp_indices <- sample(length(p), size = nsamp,
                       replace = T, prob = p/sum(p))
samp_A <- cA[samp_indices[1:nsamp]]
samp_B <- cB[samp_indices[1:nsamp]]
# add random jitter, see BDA3 p. 76
samp_A <- samp_A + runif(nsamp, A[1] - A[2], A[2] - A[1])
samp_B <- samp_B + runif(nsamp, B[1] - B[2], B[2] - B[1])

Compute LD50 conditional beta > 0

bpi <- samp_B > 0
samp_ld50 <- -samp_A[bpi]/samp_B[bpi]

Create a plot of the posterior density

# limits for the plots
xl <- c(-1.5, 7)
yl <- c(-5, 35)
pos <- ggplot(data = data.frame(cA ,cB, p), aes(x = cA, y = cB)) +
  geom_raster(aes(fill = p, alpha = p), interpolate = T) +
  geom_contour(aes(z = p), colour = 'black', size = 0.2) +
  coord_cartesian(xlim = xl, ylim = yl) +
  labs(x = 'alpha', y = 'beta') +
  scale_fill_gradient(low = 'yellow', high = 'red', guide = F) +
  scale_alpha(range = c(0, 1), guide = F)

Plot of the samples

sam <- ggplot(data = data.frame(samp_A, samp_B)) +
  geom_point(aes(samp_A, samp_B), color = 'blue', size = 0.3) +
  coord_cartesian(xlim = xl, ylim = yl) +
  labs(x = 'alpha', y = 'beta')

Plot of the histogram of LD50

his <- ggplot() +
  geom_histogram(aes(samp_ld50), binwidth = 0.04,
                 fill = 'steelblue', color = 'black') +
  coord_cartesian(xlim = c(-0.8, 0.8)) +
  labs(x = 'LD50 = -alpha/beta')

Normal approximation for Bioassay model.

Define the function to be optimized

bioassayfun <- function(w, df) {
  z <- w[1] + w[2]*df$x
  -sum(df$y*(z) - df$n*log1p(exp(z)))
}

Optimize

w0 <- c(0,0)
optim_res <- optim(w0, bioassayfun, gr = NULL, df1, hessian = T)
w <- optim_res$par
S <- solve(optim_res$hessian)

Multivariate normal probability density function

dmvnorm <- function(x, mu, sig)
  exp(-0.5*(length(x)*log(2*pi) + log(det(sig)) + (x-mu)%*%solve(sig, x-mu)))

Evaluate likelihood at points (cA,cB) this is just for illustration and would not be needed otherwise

p <- apply(cbind(cA, cB), 1, dmvnorm, w, S)

# sample from the multivariate normal 
normsamp <- mvrnorm(nsamp, w, S)

Samples of LD50 conditional beta > 0: Normal approximation does not take into account that the posterior is not symmetric and that there is very low density for negative beta values. Based on the draws from the normal approximation is is estimated that there is about 5% probability that beta is negative!

bpi <- normsamp[,2] > 0
normsamp_ld50 <- -normsamp[bpi,1]/normsamp[bpi,2]

Create a plot of the posterior density

pos_norm <- ggplot(data = data.frame(cA ,cB, p), aes(x = cA, y = cB)) +
  geom_raster(aes(fill = p, alpha = p), interpolate = T) +
  geom_contour(aes(z = p), colour = 'black', size = 0.2) +
  coord_cartesian(xlim = xl, ylim = yl) +
  labs(x = 'alpha', y = 'beta') +
  scale_fill_gradient(low = 'yellow', high = 'red', guide = F) +
  scale_alpha(range = c(0, 1), guide = F)

Plot of the samples

sam_norm <- ggplot(data = data.frame(samp_A=normsamp[,1], samp_B=normsamp[,2])) +
  geom_point(aes(samp_A, samp_B), color = 'blue', size = 0.3) +
  coord_cartesian(xlim = xl, ylim = yl) +
  labs(x = 'alpha', y = 'beta')

Plot of the histogram of LD50

his_norm <- ggplot() +
  geom_histogram(aes(normsamp_ld50), binwidth = 0.04,
                 fill = 'steelblue', color = 'black') +
  coord_cartesian(xlim = c(-0.8, 0.8)) +
  labs(x = 'LD50 = -alpha/beta, beta > 0')

Combine the plots

grid.arrange(pos, sam, his, pos_norm, sam_norm, his_norm, ncol = 3)
## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please use `guide = "none"` instead.
## It is deprecated to specify `guide = FALSE` to remove a guide. Please use `guide = "none"` instead.
## It is deprecated to specify `guide = FALSE` to remove a guide. Please use `guide = "none"` instead.
## It is deprecated to specify `guide = FALSE` to remove a guide. Please use `guide = "none"` instead.

LS0tCnRpdGxlOiAiQmF5ZXNpYW4gZGF0YSBhbmFseXNpcyBkZW1vIDQuMSIKYXV0aG9yOiAiQWtpIFZlaHRhcmksIE1hcmt1cyBQYWFzaW5pZW1pIgpkYXRlOiAiYHIgZm9ybWF0KFN5cy5EYXRlKCkpYCIKb3V0cHV0OgogIGh0bWxfZG9jdW1lbnQ6CiAgICB0aGVtZTogcmVhZGFibGUKICAgIGNvZGVfZG93bmxvYWQ6IHRydWUKLS0tCiMjIE5vcm1hbCBhcHByb3hpbWF0aW9uIGZvciBCaW9hc3NheSBtb2RlbC4KCmdncGxvdDIsIGdyaWQsIGFuZCBncmlkRXh0cmEgYXJlIHVzZWQgZm9yIHBsb3R0aW5nLCB0aWR5ciBmb3IKbWFuaXB1bGF0aW5nIGRhdGEgZnJhbWVzCgpgYGB7ciBzZXR1cCwgbWVzc2FnZT1GQUxTRSwgZXJyb3I9RkFMU0UsIHdhcm5pbmc9RkFMU0V9CmxpYnJhcnkoZ2dwbG90MikKbGlicmFyeShncmlkRXh0cmEpCmxpYnJhcnkodGlkeXIpCmxpYnJhcnkoTUFTUykKYGBgCgpCaW9hc3NheSBkYXRhLCAoQkRBMyBwYWdlIDg2KQoKYGBge3IgfQpkZjEgPC0gZGF0YS5mcmFtZSgKICB4ID0gYygtMC44NiwgLTAuMzAsIC0wLjA1LCAwLjczKSwKICBuID0gYyg1LCA1LCA1LCA1KSwKICB5ID0gYygwLCAxLCAzLCA1KQopCmBgYAoKIyMjIEdyaWQgc2FtcGxpbmcgZm9yIEJpb2Fzc2F5IG1vZGVsLgpDb21wdXRlIHRoZSBwb3N0ZXJpb3IgZGVuc2l0eSBpbiBhIGdyaWQKCi0gdXN1YWxseSBzaG91bGQgYmUgY29tcHV0ZWQgaW4gbG9nYXJpdGhtcyEKLSB3aXRoIGFsdGVybmF0aXZlIHByaW9yLCBjaGVjayB0aGF0IHJhbmdlIGFuZCBzcGFjaW5nIG9mIEEgYW5kIEIKICBhcmUgc2Vuc2libGUKCmBgYHtyIH0KQSA9IHNlcSgtMS41LCA3LCBsZW5ndGgub3V0ID0gMTAwKQpCID0gc2VxKC01LCAzNSwgbGVuZ3RoLm91dCA9IDEwMCkKIyBtYWtlIHZlY3RvcnMgdGhhdCBjb250YWluIGFsbCBwYWlyd2lzZSBjb21iaW5hdGlvbnMgb2YgQSBhbmQgQgpjQSA8LSByZXAoQSwgZWFjaCA9IGxlbmd0aChCKSkKY0IgPC0gcmVwKEIsIGxlbmd0aChBKSkKIyBhIGhlbHBlciBmdW5jdGlvbiB0byBjYWxjdWxhdGUgdGhlIGxvZyBsaWtlbGlob29kCmxvZ2wgPC0gZnVuY3Rpb24oZGYsIGEsIGIpCiAgZGZbJ3knXSooYSArIGIqZGZbJ3gnXSkgLSBkZlsnbiddKmxvZzFwKGV4cChhICsgYipkZlsneCddKSkKIyBjYWxjdWxhdGUgbGlrZWxpaG9vZHM6IGFwcGx5IGxvZ2wgZnVuY3Rpb24gZm9yIGVhY2ggb2JzZXJ2YXRpb24KIyBpZS4gZWFjaCByb3cgb2YgZGF0YSBmcmFtZSBvZiB4LCBuIGFuZCB5CnAgPC0gYXBwbHkoZGYxLCAxLCBsb2dsLCBjQSwgY0IpICU+JSByb3dTdW1zKCkgJT4lIGV4cCgpCmBgYAoKU2FtcGxlIGZyb20gdGhlIGdyaWQgKHdpdGggcmVwbGFjZW1lbnQpCgpgYGB7ciB9Cm5zYW1wIDwtIDEwMDAKc2FtcF9pbmRpY2VzIDwtIHNhbXBsZShsZW5ndGgocCksIHNpemUgPSBuc2FtcCwKICAgICAgICAgICAgICAgICAgICAgICByZXBsYWNlID0gVCwgcHJvYiA9IHAvc3VtKHApKQpzYW1wX0EgPC0gY0Fbc2FtcF9pbmRpY2VzWzE6bnNhbXBdXQpzYW1wX0IgPC0gY0Jbc2FtcF9pbmRpY2VzWzE6bnNhbXBdXQojIGFkZCByYW5kb20gaml0dGVyLCBzZWUgQkRBMyBwLiA3NgpzYW1wX0EgPC0gc2FtcF9BICsgcnVuaWYobnNhbXAsIEFbMV0gLSBBWzJdLCBBWzJdIC0gQVsxXSkKc2FtcF9CIDwtIHNhbXBfQiArIHJ1bmlmKG5zYW1wLCBCWzFdIC0gQlsyXSwgQlsyXSAtIEJbMV0pCmBgYAoKQ29tcHV0ZSBMRDUwIGNvbmRpdGlvbmFsIGJldGEgPiAwCgpgYGB7ciB9CmJwaSA8LSBzYW1wX0IgPiAwCnNhbXBfbGQ1MCA8LSAtc2FtcF9BW2JwaV0vc2FtcF9CW2JwaV0KYGBgCgpDcmVhdGUgYSBwbG90IG9mIHRoZSBwb3N0ZXJpb3IgZGVuc2l0eQoKYGBge3IgfQojIGxpbWl0cyBmb3IgdGhlIHBsb3RzCnhsIDwtIGMoLTEuNSwgNykKeWwgPC0gYygtNSwgMzUpCnBvcyA8LSBnZ3Bsb3QoZGF0YSA9IGRhdGEuZnJhbWUoY0EgLGNCLCBwKSwgYWVzKHggPSBjQSwgeSA9IGNCKSkgKwogIGdlb21fcmFzdGVyKGFlcyhmaWxsID0gcCwgYWxwaGEgPSBwKSwgaW50ZXJwb2xhdGUgPSBUKSArCiAgZ2VvbV9jb250b3VyKGFlcyh6ID0gcCksIGNvbG91ciA9ICdibGFjaycsIHNpemUgPSAwLjIpICsKICBjb29yZF9jYXJ0ZXNpYW4oeGxpbSA9IHhsLCB5bGltID0geWwpICsKICBsYWJzKHggPSAnYWxwaGEnLCB5ID0gJ2JldGEnKSArCiAgc2NhbGVfZmlsbF9ncmFkaWVudChsb3cgPSAneWVsbG93JywgaGlnaCA9ICdyZWQnLCBndWlkZSA9IEYpICsKICBzY2FsZV9hbHBoYShyYW5nZSA9IGMoMCwgMSksIGd1aWRlID0gRikKYGBgCgpQbG90IG9mIHRoZSBzYW1wbGVzCgpgYGB7ciB9CnNhbSA8LSBnZ3Bsb3QoZGF0YSA9IGRhdGEuZnJhbWUoc2FtcF9BLCBzYW1wX0IpKSArCiAgZ2VvbV9wb2ludChhZXMoc2FtcF9BLCBzYW1wX0IpLCBjb2xvciA9ICdibHVlJywgc2l6ZSA9IDAuMykgKwogIGNvb3JkX2NhcnRlc2lhbih4bGltID0geGwsIHlsaW0gPSB5bCkgKwogIGxhYnMoeCA9ICdhbHBoYScsIHkgPSAnYmV0YScpCmBgYAoKUGxvdCBvZiB0aGUgaGlzdG9ncmFtIG9mIExENTAKCmBgYHtyIH0KaGlzIDwtIGdncGxvdCgpICsKICBnZW9tX2hpc3RvZ3JhbShhZXMoc2FtcF9sZDUwKSwgYmlud2lkdGggPSAwLjA0LAogICAgICAgICAgICAgICAgIGZpbGwgPSAnc3RlZWxibHVlJywgY29sb3IgPSAnYmxhY2snKSArCiAgY29vcmRfY2FydGVzaWFuKHhsaW0gPSBjKC0wLjgsIDAuOCkpICsKICBsYWJzKHggPSAnTEQ1MCA9IC1hbHBoYS9iZXRhJykKYGBgCgojIyMgTm9ybWFsIGFwcHJveGltYXRpb24gZm9yIEJpb2Fzc2F5IG1vZGVsLgpEZWZpbmUgdGhlIGZ1bmN0aW9uIHRvIGJlIG9wdGltaXplZAoKYGBge3IgfQpiaW9hc3NheWZ1biA8LSBmdW5jdGlvbih3LCBkZikgewogIHogPC0gd1sxXSArIHdbMl0qZGYkeAogIC1zdW0oZGYkeSooeikgLSBkZiRuKmxvZzFwKGV4cCh6KSkpCn0KYGBgCgpPcHRpbWl6ZQoKYGBge3IgfQp3MCA8LSBjKDAsMCkKb3B0aW1fcmVzIDwtIG9wdGltKHcwLCBiaW9hc3NheWZ1biwgZ3IgPSBOVUxMLCBkZjEsIGhlc3NpYW4gPSBUKQp3IDwtIG9wdGltX3JlcyRwYXIKUyA8LSBzb2x2ZShvcHRpbV9yZXMkaGVzc2lhbikKYGBgCgpNdWx0aXZhcmlhdGUgbm9ybWFsIHByb2JhYmlsaXR5IGRlbnNpdHkgZnVuY3Rpb24KCmBgYHtyIH0KZG12bm9ybSA8LSBmdW5jdGlvbih4LCBtdSwgc2lnKQogIGV4cCgtMC41KihsZW5ndGgoeCkqbG9nKDIqcGkpICsgbG9nKGRldChzaWcpKSArICh4LW11KSUqJXNvbHZlKHNpZywgeC1tdSkpKQpgYGAKCkV2YWx1YXRlIGxpa2VsaWhvb2QgYXQgcG9pbnRzIChjQSxjQikgCnRoaXMgaXMganVzdCBmb3IgaWxsdXN0cmF0aW9uIGFuZCB3b3VsZCBub3QgYmUgbmVlZGVkIG90aGVyd2lzZQoKYGBge3IgfQpwIDwtIGFwcGx5KGNiaW5kKGNBLCBjQiksIDEsIGRtdm5vcm0sIHcsIFMpCgojIHNhbXBsZSBmcm9tIHRoZSBtdWx0aXZhcmlhdGUgbm9ybWFsIApub3Jtc2FtcCA8LSBtdnJub3JtKG5zYW1wLCB3LCBTKQpgYGAKClNhbXBsZXMgb2YgTEQ1MCBjb25kaXRpb25hbCBiZXRhID4gMDoKTm9ybWFsIGFwcHJveGltYXRpb24gZG9lcyBub3QgdGFrZSBpbnRvIGFjY291bnQgdGhhdCB0aGUgcG9zdGVyaW9yCmlzIG5vdCBzeW1tZXRyaWMgYW5kIHRoYXQgdGhlcmUgaXMgdmVyeSBsb3cgZGVuc2l0eSBmb3IgbmVnYXRpdmUKYmV0YSB2YWx1ZXMuIEJhc2VkIG9uIHRoZSBkcmF3cyBmcm9tIHRoZSBub3JtYWwgYXBwcm94aW1hdGlvbgppcyBpcyBlc3RpbWF0ZWQgdGhhdCB0aGVyZSBpcyBhYm91dCA1JSBwcm9iYWJpbGl0eSB0aGF0IGJldGEgaXMgbmVnYXRpdmUhCgpgYGB7ciB9CmJwaSA8LSBub3Jtc2FtcFssMl0gPiAwCm5vcm1zYW1wX2xkNTAgPC0gLW5vcm1zYW1wW2JwaSwxXS9ub3Jtc2FtcFticGksMl0KYGBgCgpDcmVhdGUgYSBwbG90IG9mIHRoZSBwb3N0ZXJpb3IgZGVuc2l0eQoKYGBge3IgfQpwb3Nfbm9ybSA8LSBnZ3Bsb3QoZGF0YSA9IGRhdGEuZnJhbWUoY0EgLGNCLCBwKSwgYWVzKHggPSBjQSwgeSA9IGNCKSkgKwogIGdlb21fcmFzdGVyKGFlcyhmaWxsID0gcCwgYWxwaGEgPSBwKSwgaW50ZXJwb2xhdGUgPSBUKSArCiAgZ2VvbV9jb250b3VyKGFlcyh6ID0gcCksIGNvbG91ciA9ICdibGFjaycsIHNpemUgPSAwLjIpICsKICBjb29yZF9jYXJ0ZXNpYW4oeGxpbSA9IHhsLCB5bGltID0geWwpICsKICBsYWJzKHggPSAnYWxwaGEnLCB5ID0gJ2JldGEnKSArCiAgc2NhbGVfZmlsbF9ncmFkaWVudChsb3cgPSAneWVsbG93JywgaGlnaCA9ICdyZWQnLCBndWlkZSA9IEYpICsKICBzY2FsZV9hbHBoYShyYW5nZSA9IGMoMCwgMSksIGd1aWRlID0gRikKYGBgCgpQbG90IG9mIHRoZSBzYW1wbGVzCgpgYGB7ciB9CnNhbV9ub3JtIDwtIGdncGxvdChkYXRhID0gZGF0YS5mcmFtZShzYW1wX0E9bm9ybXNhbXBbLDFdLCBzYW1wX0I9bm9ybXNhbXBbLDJdKSkgKwogIGdlb21fcG9pbnQoYWVzKHNhbXBfQSwgc2FtcF9CKSwgY29sb3IgPSAnYmx1ZScsIHNpemUgPSAwLjMpICsKICBjb29yZF9jYXJ0ZXNpYW4oeGxpbSA9IHhsLCB5bGltID0geWwpICsKICBsYWJzKHggPSAnYWxwaGEnLCB5ID0gJ2JldGEnKQpgYGAKClBsb3Qgb2YgdGhlIGhpc3RvZ3JhbSBvZiBMRDUwCgpgYGB7ciB9Cmhpc19ub3JtIDwtIGdncGxvdCgpICsKICBnZW9tX2hpc3RvZ3JhbShhZXMobm9ybXNhbXBfbGQ1MCksIGJpbndpZHRoID0gMC4wNCwKICAgICAgICAgICAgICAgICBmaWxsID0gJ3N0ZWVsYmx1ZScsIGNvbG9yID0gJ2JsYWNrJykgKwogIGNvb3JkX2NhcnRlc2lhbih4bGltID0gYygtMC44LCAwLjgpKSArCiAgbGFicyh4ID0gJ0xENTAgPSAtYWxwaGEvYmV0YSwgYmV0YSA+IDAnKQpgYGAKCkNvbWJpbmUgdGhlIHBsb3RzCgpgYGB7ciB9CmdyaWQuYXJyYW5nZShwb3MsIHNhbSwgaGlzLCBwb3Nfbm9ybSwgc2FtX25vcm0sIGhpc19ub3JtLCBuY29sID0gMykKYGBgCgo=