We want to get a sense of the shape of a single term's contribution to the log posterior function.
[ \log \pi(\sigma^2_e, \sigma^2_s | y) = B - \frac{1}{2}\sum_j \bigg[ c_j \log(a_j \sigma^2_s + b_j \sigma^2_e) + \frac{d_j}{(a_j \sigma^2_s + b_j \sigma^2_e)} \bigg] ]
To get a sense of the shape that each term takes, we fit a simple random effects model.
library(knitr) opts_chunk$set(echo = FALSE, warning = FALSE, message = FALSE)
library(ggplot2); library(dplyr); library(devtools); devtools::install_github("andrewpbray/lmmoptim") library(lmmoptim) data(hmo) data(hmo_states) hmobig <- left_join(hmo, hmo_states, by = "state") y <- hmobig$indPrem n <- length(y) mod <- lm(indPrem ~ expPerAdm + (region=="NE") + state, data = hmobig, x = TRUE ) x <- mod$x[, 1:3] z <- mod$x[, -(1:3), drop = FALSE] SigE <- diag(n) SigS <- diag(44) lines <- findlines(x, z, y, SigE, SigS)
In the figure below to the left, we plot only the first term in the sum, with the maximum value achieved along the dotted gray line. The value of function is mapped to color, with lowest values being red and high values being light yellow. The function steeply decends near the origin, so that area is not color mapped so as to not saturate the color palette. The value of that the function evaluates to along the solid blue line is plotted separately below to the right, where it is a function of the linear term, $(a_j \sigma^2_s + b_j \sigma^2_e)$.
plotTerm <- function(lines, which.line = 1, npix = 100) { c_j <- lines[which.line, "multiplier.log"] d_j <- lines[which.line, "multiplier.inv"] x <- seq(0, lines[which.line, "int.sigsqe"], length.out = npix) y <- seq(0, lines[which.line, "int.sigsqs"], length.out = npix) fterm <- function(x, y, lines, which.line) { linearterm <- lines[which.line, "a"] * y + lines[which.line, "b"] * x ifelse(linearterm == 0, NA, -0.5 * (lines[which.line, "multiplier.log"] * log(linearterm) + lines[which.line, "multiplier.inv"]/linearterm)) } z <- outer(x, y, FUN = fterm, lines = lines, which.line = which.line) par(mfrow = c(1, 2)) # plot surface in original coordinates image(x, y, z, zlim = c(.05 * min(as.numeric(z), na.rm = T), max(as.numeric(z), na.rm = T))) abline(lines[which.line, "int.sigsqs"], -lines[which.line, "int.sigsqs"]/lines[which.line, "int.sigsqe"], col = "gray", lwd = 2, lty = 2) perpslope <- tan(pi/2 - atan(lines[which.line, "int.sigsqs"]/lines[which.line, "int.sigsqe"])) abline(0, 1/perpslope, col = "steelblue", lwd = 2) # plot trace of linear component linearterm <- lines[which.line, "a"] * y + x plot(linearterm, -0.5 * (c_j * log(linearterm) + d_j/linearterm), type = "n", ylab = "f", xlab = "linear term") lines(linearterm, -0.5 * (c_j * log(linearterm) + d_j/linearterm), col = "steelblue", lwd = 2) abline(v = linearterm[which.max(-0.5 * (c_j * log(linearterm) + d_j/linearterm))], col = "tomato") } plotTerm(lines, which.line = 1, npix = 200)
We can see that the function drops off very sharply towards the origin and very gradually away from it. It's maximum is indicated with the vertical red line.
Each term in the sum will have it's own set of line parameters, ${ a_j, b_j}$, and function shape parameters, ${c_j, d_j }$. $c_j$ is one for all terms except that corresponding to the prior, while the values of $d_j$ range the first three orders of magnitude (at least in the HMO model). Below we plot the shape of the function corresponding to three different values of $d_j$ and indicate the max with a vertical red line (note the differing y axes).
x <- 1:6000 minv <- c(30, 300, 3000) maxes <- rep(NA, 3) par(mfrow = c(1, 3)) for(i in 1:3) { fx <- -0.5 * (log(x) + minv[i]/x) plot(x, fx, type = "n", ylab = "f", xlab = "linear term", main = bquote(d[j] ~ "=" ~ .(minv[i]))) lines(x, fx, col = "steelblue", lwd = 2) abline(v = which.max(fx), col = "tomato") maxes[i] <- max(fx) }
Some take away lessons from this:
r round(maxes, 2)
.Can we be more efficient by addressing the terms in order by descending $d_j$? Is there any sensible way to rescale each term so that some of these parameters are standardized across the terms, then optimize, then undo the scaling?
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.