Safeguarded Newton-Raphson for Logistic Regression
#| '!! shinylive warning !!': |
#| shinylive does not work in self-contained HTML documents.
#| Please set `embed-resources: false` in your metadata.
#| standalone: true
#| viewerHeight: 960
library(shiny)
### ---------- fixed data ----------------------------------------------------
set.seed(1)
n <- 200
z <- sort(runif(n, -2, 2))
r <- rbinom(n, 1, plogis(0 + 1.5 * z))
r_jit <- r + runif(n, -0.1, 0.1) # vertical jitter, for display only
log1pexp <- function(u) pmax(u, 0) + log1p(exp(-abs(u)))
negll <- function(b) {
u <- b[1] + b[2] * z
-sum(u * r - log1pexp(u))
}
## gradient of the negative log-likelihood (minus the score)
grad_negll <- function(b) {
p <- plogis(b[1] + b[2] * z)
-c(sum(r - p), sum(z * (r - p)))
}
## observed information = Hessian of the negative log-likelihood
info_mat <- function(b) {
p <- plogis(b[1] + b[2] * z)
w <- p * (1 - p)
matrix(c(sum(w), sum(z * w), sum(z * w), sum(z * z * w)), 2, 2)
}
cond_num <- function(J) {
e <- abs(eigen(J, symmetric = TRUE, only.values = TRUE)$values)
if (min(e) <= 0) Inf else max(e) / min(e)
}
### ---------- robustification of the Hessian --------------------------------
## Replace J by a nearby positive definite matrix whose condition number is at
## most 1/tau, so that J^{-1} stays bounded and the step remains a descent
## direction. "ridge" adds lambda*I (Levenberg-Marquardt); "eigen" raises only
## the eigenvalues that are too small (eigenvalue flooring).
robustify <- function(J, method, tau = 1e-3) {
if (method == "none" || !all(is.finite(J))) return(J)
ev <- eigen(J, symmetric = TRUE)
lmax <- max(ev$values)
if (!is.finite(lmax) || lmax <= 0) return(diag(2))
if (method == "ridge") {
J + max(0, tau * lmax - min(ev$values)) * diag(2)
} else {
ev$vectors %*% (pmax(ev$values, tau * lmax) * t(ev$vectors))
}
}
### ---------- contour grid (computed once) ----------------------------------
B0 <- seq(-3, 3, by = 0.1)
B1 <- seq(-8, 8, by = 0.2)
Z <- outer(B0, B1, Vectorize(function(a, b) negll(c(a, b))))
### ---------- Newton-Raphson path -------------------------------------------
## Returns the iterates with their negll, gradient, accepted step length and the
## condition number of the (unmodified) information matrix, stopping early on
## convergence or when the iteration breaks down.
nr_path <- function(b0, maxit = 15, ls = "none", robust = "none", tol = 1e-8) {
cn <- c("beta0", "beta1", "negll", "g0", "g1", "s", "kappa")
P <- matrix(NA_real_, maxit + 1, length(cn), dimnames = list(NULL, cn))
b <- b0
P[1, ] <- c(b, negll(b), grad_negll(b), NA, cond_num(info_mat(b)))
n_ok <- 1; note <- "maximum number of iterations reached"
for (i in seq_len(maxit)) {
g <- P[i, 4:5]
J <- robustify(info_mat(b), robust)
step <- tryCatch(solve(J, g), error = function(e) rep(NA_real_, 2))
if (any(!is.finite(step))) {
note <- "information matrix is numerically singular: stopped"
break
}
## backtracking line search: halve the step until the Armijo condition
## negll(b - s*step) <= negll(b) - c1*s*<g, step> is satisfied
s <- 1
if (ls == "armijo") {
f0 <- P[i, 3]; dd <- sum(g * step)
while (s > 1e-10) {
fn <- negll(b - s * step)
if (is.finite(fn) && fn <= f0 - 1e-4 * s * dd) break
s <- s / 2
}
if (s <= 1e-10) {
note <- "line search found no decrease: stopped"
break
}
}
b <- b - s * step
if (any(!is.finite(b))) { note <- "iterate is no longer finite: diverged"; break }
P[i + 1, ] <- c(b, negll(b), grad_negll(b), s, cond_num(info_mat(b)))
n_ok <- i + 1
if (max(abs(s * step)) < tol) { note <- "converged"; break }
}
list(P = P[seq_len(n_ok), , drop = FALSE], n = n_ok, note = note)
}
### ---------- UI -------------------------------------------------------------
## an icon-only playback button; the title shows as a tooltip
ctrl_btn <- function(id, icon_name, title, class = "btn-default")
actionButton(id, NULL, icon = icon(icon_name), title = title, class = class,
style = "flex: 1; padding: 6px 0;")
ui <- fluidPage(
## control panel
wellPanel(
fluidRow(
column(4, fluidRow(
column(6, numericInput("b0", HTML("starting β<sub>0</sub>"),
value = 0, step = 0.5, width = "100%")),
column(6, numericInput("b1", HTML("starting β<sub>1</sub>"),
value = 5, step = 0.5, width = "100%"))),
helpText("Or click the contour plot to choose a new start.",
style = "margin-top: -8px;")),
column(4, selectInput("ls", "line search",
c("none (full Newton step)" = "none",
"backtracking (Armijo)" = "armijo"),
width = "100%")),
column(4, selectInput("robust", "robustification of the Hessian",
c("none (observed information)" = "none",
"ridge (Levenberg-Marquardt)" = "ridge",
"eigenvalue flooring" = "eigen"),
width = "100%"))
),
fluidRow(
column(4, div(style = "display: flex; gap: 4px; padding-top: 25px;",
ctrl_btn("first", "backward-fast", "Back to the start"),
ctrl_btn("back", "backward-step", "One step back"),
ctrl_btn("play", "play", "Play", class = "btn-primary"),
ctrl_btn("pause", "pause", "Pause"),
ctrl_btn("step", "forward-step", "One step forward"),
ctrl_btn("toend", "forward-fast", "To the end"))),
column(3, div(style = "padding-top: 12px;",
checkboxInput("compare", "overlay plain Newton-Raphson",
TRUE),
checkboxInput("autoplay", "autoplay on a new start",
TRUE))),
column(5, sliderInput("iter", "iteration", min = 0, max = 1,
value = 0, step = 1, width = "100%"))
)
),
## the two views of the same iterate
fluidRow(
column(6, plotOutput("contour", click = "click", height = "430px")),
column(6, plotOutput("fit", height = "430px"))
),
## details of the iterations, shown under the plots
hr(),
fluidRow(
column(7, tableOutput("tab"), htmlOutput("note")),
column(5, helpText("Try the starting values (0, 3) and then (0, 5), first ",
"with both safeguards off, then with each in turn. From ",
"(0, 5) the plain iteration stops because the ",
"information matrix goes numerically singular; the ",
"robustification is what keeps it invertible."))
)
)
### ---------- server ---------------------------------------------------------
server <- function(input, output, session) {
playing <- reactiveVal(FALSE)
play_ms <- 900
## clicking the contour plot moves the starting value
observeEvent(input$click, {
p <- c(input$click$x, input$click$y)
if (all(is.finite(p))) {
updateNumericInput(session, "b0", value = round(p[1], 2))
updateNumericInput(session, "b1", value = round(p[2], 2))
}
})
start <- reactive({
req(is.finite(input$b0), is.finite(input$b1))
c(input$b0, input$b1)
})
modified <- reactive(input$ls != "none" || input$robust != "none")
path <- reactive(nr_path(start(), maxit = 15, ls = input$ls,
robust = input$robust))
## the plain (unsafeguarded) path, shown for comparison
ref <- reactive({
if (!isTRUE(input$compare) || !modified()) return(NULL)
nr_path(start(), maxit = 15)
})
## current iteration, clamped to the length of the path
k <- reactive(min(as.integer(input$iter), path()$n - 1))
## a new path restarts the display at iteration 0
observeEvent(path(), {
updateSliderInput(session, "iter", max = max(1, path()$n - 1), value = 0)
playing(isTRUE(input$autoplay) && path()$n > 1)
})
go_to <- function(i) {
updateSliderInput(session, "iter",
value = max(0, min(i, path()$n - 1)))
}
observeEvent(input$play, {
if (k() >= path()$n - 1) go_to(0)
if (path()$n > 1) playing(TRUE)
})
observeEvent(input$pause, playing(FALSE))
observeEvent(input$step, { playing(FALSE); go_to(k() + 1) })
observeEvent(input$back, { playing(FALSE); go_to(k() - 1) })
observeEvent(input$first, { playing(FALSE); go_to(0) })
observeEvent(input$toend, { playing(FALSE); go_to(path()$n - 1) })
observe({
if (!isTRUE(playing())) return()
kk <- isolate(k())
if (kk >= isolate(path()$n) - 1) { playing(FALSE); return() }
invalidateLater(play_ms, session)
isolate(go_to(kk + 1))
})
output$contour <- renderPlot({
P <- path()$P; kk <- k() + 1L
RP <- ref()
par(mar = c(4, 4, 2, 1))
contour(B0, B1, Z, nlevels = 40, drawlabels = FALSE,
col = "grey70",
xlab = expression(beta[0]), ylab = expression(beta[1]),
main = "negative log-likelihood surface and NR path")
points(0, 1.5, pch = 3, col = "darkgreen", lwd = 2)
## the unsafeguarded path, drawn underneath for comparison
if (!is.null(RP)) {
Q <- RP$P[seq_len(min(kk, RP$n)), , drop = FALSE]
if (nrow(Q) > 1) lines(Q[, 1], Q[, 2], col = "grey45", lwd = 2,
lty = 2)
points(Q[, 1], Q[, 2], pch = 21, bg = "grey85", col = "grey45",
cex = 0.9)
}
if (kk > 1) lines(P[1:kk, 1], P[1:kk, 2], col = "steelblue", lwd = 2)
points(P[1:kk, 1], P[1:kk, 2], pch = 21, bg = "steelblue", cex = 1.1)
cur <- P[kk, ]
inside <- cur[1] >= min(B0) && cur[1] <= max(B0) &&
cur[2] >= min(B1) && cur[2] <= max(B1)
if (inside) {
## Arrow along the score (uphill in the log-likelihood, i.e. minus
## the gradient of negll), of fixed length on screen. The gradient
## is a covector: because the two axes are drawn at different
## scales, its on-screen components are obtained by MULTIPLYING by
## the units-per-inch factors (not dividing, as for a displacement).
## Only then is the arrow perpendicular to the contours as drawn.
g <- -cur[4:5]
upi <- c(diff(par("usr")[1:2]) / par("pin")[1],
diff(par("usr")[3:4]) / par("pin")[2]) # units per inch
gs <- g * upi # inches
## at a stationary point the direction is numerical noise, so the
## arrow is dropped rather than drawn at full length
if (sqrt(sum(g^2)) > 1e-6 && sqrt(sum(gs^2)) > 0) {
d <- 0.45 * gs / sqrt(sum(gs^2)) * upi # 0.45 inch long
arrows(cur[1], cur[2], cur[1] + d[1], cur[2] + d[2],
col = "#E8820C", lwd = 3, length = 0.09)
}
points(cur[1], cur[2], pch = 21, bg = "firebrick", cex = 1.8)
} else {
mtext(sprintf("current iterate (%.3g, %.3g) is off the plotted region",
cur[1], cur[2]), side = 3, line = -1.5,
col = "firebrick", cex = 0.9)
}
leg <- list(txt = c("current iterate", "true (0, 1.5)",
"gradient of log-likelihood (direction)"),
pch = c(21, 3, NA), bg = c("firebrick", NA, NA),
lwd = c(NA, 2, 3), lty = c(NA, NA, 1),
col = c("black", "darkgreen", "#E8820C"))
if (!is.null(RP)) {
leg$txt <- c(leg$txt, "plain Newton-Raphson")
leg$pch <- c(leg$pch, NA); leg$bg <- c(leg$bg, NA)
leg$lwd <- c(leg$lwd, 2); leg$lty <- c(leg$lty, 2)
leg$col <- c(leg$col, "grey45")
}
legend("bottomright", bg = "#ffffffcc", box.col = NA, cex = 0.85,
pch = leg$pch, pt.bg = leg$bg, lwd = leg$lwd, lty = leg$lty,
col = leg$col, legend = leg$txt)
})
## the logistic curve fitted at the current iterate, over the data; the 0/1
## responses are jittered vertically so that their density can be seen
output$fit <- renderPlot({
P <- path()$P; kk <- k() + 1L; RP <- ref()
xs <- seq(min(z), max(z), length.out = 200)
curve_at <- function(b) plogis(b[1] + b[2] * xs)
par(mar = c(4, 4, 2, 1))
plot(z, r_jit, pch = 19, cex = 0.7, col = adjustcolor("grey30", 0.4),
ylim = c(-0.2, 1.2), xlab = "x", ylab = "P(y = 1)",
main = sprintf("fitted logistic curve, iteration %d", kk - 1))
abline(h = c(0, 1), col = "grey85")
lines(xs, curve_at(c(0, 1.5)), col = "darkgreen", lwd = 2, lty = 3)
## the earlier iterates, faint
if (kk > 1) for (i in seq_len(kk - 1))
lines(xs, curve_at(P[i, 1:2]), col = adjustcolor("steelblue", 0.35))
## the plain (unsafeguarded) iterate at the same step
if (!is.null(RP))
lines(xs, curve_at(RP$P[min(kk, RP$n), 1:2]), col = "grey45",
lwd = 2, lty = 2)
lines(xs, curve_at(P[kk, 1:2]), col = "firebrick", lwd = 3)
leg <- list(txt = c("data (jittered)", "current fit", "earlier iterates",
"true curve"),
pch = c(19, NA, NA, NA), lwd = c(NA, 3, 1, 2),
lty = c(NA, 1, 1, 3),
col = c(adjustcolor("grey30", 0.6), "firebrick", "steelblue",
"darkgreen"))
if (!is.null(RP)) {
leg$txt <- c(leg$txt, "plain Newton-Raphson")
leg$pch <- c(leg$pch, NA); leg$lwd <- c(leg$lwd, 2)
leg$lty <- c(leg$lty, 2); leg$col <- c(leg$col, "grey45")
}
legend("right", bg = "#ffffffcc", box.col = NA, cex = 0.85,
pch = leg$pch, lwd = leg$lwd, lty = leg$lty, col = leg$col,
legend = leg$txt)
})
## one decimal throughout; a diverging run reaches values that would
## otherwise be far too wide, so those switch to one-decimal scientific
f1 <- function(x) ifelse(!is.finite(x), "",
ifelse(abs(x) >= 1e5 | (x != 0 & abs(x) < 1e-3),
formatC(x, format = "e", digits = 1),
formatC(x, format = "f", digits = 1)))
## why the run ended, shown as one line under the table
output$note <- renderUI({
R <- path(); RP <- ref()
if (k() + 1 < R$n) return(NULL)
HTML(paste0("<i>", R$note,
if (!is.null(RP)) paste0("; plain Newton-Raphson: ",
RP$note),
"</i>"))
})
output$tab <- renderTable({
P <- path()$P[seq_len(k() + 1), , drop = FALSE]
d <- data.frame(iteration = as.character(seq_len(nrow(P)) - 1),
beta0 = f1(P[, 1]), beta1 = f1(P[, 2]),
negll = f1(P[, 3]), step = f1(P[, 6]),
`cond(J)` = f1(P[, 7]),
check.names = FALSE, stringsAsFactors = FALSE)
tail(d, 5)
}, rownames = FALSE, align = "r", width = "100%")
}
shinyApp(ui, server)About the app
Shows how the full Newton-Raphson step diverges on a logistic-regression likelihood started far from the maximum, and how backtracking on the step length rescues the same starting value. The iterates are traced over the contours of the negative log-likelihood, with the orange arrow marking the gradient direction at the current iterate and the unsafeguarded path left in place as a dashed grey curve for comparison.
Click in the contour plot (or type values) to choose a starting point, then use the playback buttons to step through the iterations or play the path as an animation over the contours of the negative log-likelihood. The orange arrow shows the direction of the gradient of the log-likelihood. Switching on either safeguard keeps the plain iteration visible as a dashed grey path for comparison; the details of each iteration, and the reason the run ended, are reported under the plots. The right-hand plot shows the logistic curve fitted at the current iterate over the data, whose 0/1 responses are jittered vertically so that their density can be seen.
NoteR source for this app
library(shiny)
### ---------- fixed data ----------------------------------------------------
set.seed(1)
n <- 200
z <- sort(runif(n, -2, 2))
r <- rbinom(n, 1, plogis(0 + 1.5 * z))
r_jit <- r + runif(n, -0.1, 0.1) # vertical jitter, for display only
log1pexp <- function(u) pmax(u, 0) + log1p(exp(-abs(u)))
negll <- function(b) {
u <- b[1] + b[2] * z
-sum(u * r - log1pexp(u))
}
## gradient of the negative log-likelihood (minus the score)
grad_negll <- function(b) {
p <- plogis(b[1] + b[2] * z)
-c(sum(r - p), sum(z * (r - p)))
}
## observed information = Hessian of the negative log-likelihood
info_mat <- function(b) {
p <- plogis(b[1] + b[2] * z)
w <- p * (1 - p)
matrix(c(sum(w), sum(z * w), sum(z * w), sum(z * z * w)), 2, 2)
}
cond_num <- function(J) {
e <- abs(eigen(J, symmetric = TRUE, only.values = TRUE)$values)
if (min(e) <= 0) Inf else max(e) / min(e)
}
### ---------- robustification of the Hessian --------------------------------
## Replace J by a nearby positive definite matrix whose condition number is at
## most 1/tau, so that J^{-1} stays bounded and the step remains a descent
## direction. "ridge" adds lambda*I (Levenberg-Marquardt); "eigen" raises only
## the eigenvalues that are too small (eigenvalue flooring).
robustify <- function(J, method, tau = 1e-3) {
if (method == "none" || !all(is.finite(J))) return(J)
ev <- eigen(J, symmetric = TRUE)
lmax <- max(ev$values)
if (!is.finite(lmax) || lmax <= 0) return(diag(2))
if (method == "ridge") {
J + max(0, tau * lmax - min(ev$values)) * diag(2)
} else {
ev$vectors %*% (pmax(ev$values, tau * lmax) * t(ev$vectors))
}
}
### ---------- contour grid (computed once) ----------------------------------
B0 <- seq(-3, 3, by = 0.1)
B1 <- seq(-8, 8, by = 0.2)
Z <- outer(B0, B1, Vectorize(function(a, b) negll(c(a, b))))
### ---------- Newton-Raphson path -------------------------------------------
## Returns the iterates with their negll, gradient, accepted step length and the
## condition number of the (unmodified) information matrix, stopping early on
## convergence or when the iteration breaks down.
nr_path <- function(b0, maxit = 15, ls = "none", robust = "none", tol = 1e-8) {
cn <- c("beta0", "beta1", "negll", "g0", "g1", "s", "kappa")
P <- matrix(NA_real_, maxit + 1, length(cn), dimnames = list(NULL, cn))
b <- b0
P[1, ] <- c(b, negll(b), grad_negll(b), NA, cond_num(info_mat(b)))
n_ok <- 1; note <- "maximum number of iterations reached"
for (i in seq_len(maxit)) {
g <- P[i, 4:5]
J <- robustify(info_mat(b), robust)
step <- tryCatch(solve(J, g), error = function(e) rep(NA_real_, 2))
if (any(!is.finite(step))) {
note <- "information matrix is numerically singular: stopped"
break
}
## backtracking line search: halve the step until the Armijo condition
## negll(b - s*step) <= negll(b) - c1*s*<g, step> is satisfied
s <- 1
if (ls == "armijo") {
f0 <- P[i, 3]; dd <- sum(g * step)
while (s > 1e-10) {
fn <- negll(b - s * step)
if (is.finite(fn) && fn <= f0 - 1e-4 * s * dd) break
s <- s / 2
}
if (s <= 1e-10) {
note <- "line search found no decrease: stopped"
break
}
}
b <- b - s * step
if (any(!is.finite(b))) { note <- "iterate is no longer finite: diverged"; break }
P[i + 1, ] <- c(b, negll(b), grad_negll(b), s, cond_num(info_mat(b)))
n_ok <- i + 1
if (max(abs(s * step)) < tol) { note <- "converged"; break }
}
list(P = P[seq_len(n_ok), , drop = FALSE], n = n_ok, note = note)
}
### ---------- UI -------------------------------------------------------------
## an icon-only playback button; the title shows as a tooltip
ctrl_btn <- function(id, icon_name, title, class = "btn-default")
actionButton(id, NULL, icon = icon(icon_name), title = title, class = class,
style = "flex: 1; padding: 6px 0;")
ui <- fluidPage(
titlePanel("Shinylive App for Safeguarded Newton-Raphson"),
## control panel
wellPanel(
fluidRow(
column(4, fluidRow(
column(6, numericInput("b0", HTML("starting β<sub>0</sub>"),
value = 0, step = 0.5, width = "100%")),
column(6, numericInput("b1", HTML("starting β<sub>1</sub>"),
value = 5, step = 0.5, width = "100%"))),
helpText("Or click the contour plot to choose a new start.",
style = "margin-top: -8px;")),
column(4, selectInput("ls", "line search",
c("none (full Newton step)" = "none",
"backtracking (Armijo)" = "armijo"),
width = "100%")),
column(4, selectInput("robust", "robustification of the Hessian",
c("none (observed information)" = "none",
"ridge (Levenberg-Marquardt)" = "ridge",
"eigenvalue flooring" = "eigen"),
width = "100%"))
),
fluidRow(
column(5, sliderInput("iter", "iteration", min = 0, max = 1,
value = 0, step = 1, width = "100%")),
column(4, div(style = "display: flex; gap: 4px; padding-top: 25px;",
ctrl_btn("first", "backward-fast", "Back to the start"),
ctrl_btn("back", "backward-step", "One step back"),
ctrl_btn("play", "play", "Play", class = "btn-primary"),
ctrl_btn("pause", "pause", "Pause"),
ctrl_btn("step", "forward-step", "One step forward"),
ctrl_btn("toend", "forward-fast", "To the end"))),
column(3, div(style = "padding-top: 12px;",
checkboxInput("compare", "overlay plain Newton-Raphson",
TRUE),
checkboxInput("autoplay", "autoplay on a new start",
TRUE)))
)
),
## the two views of the same iterate
fluidRow(
column(6, plotOutput("contour", click = "click", height = "430px")),
column(6, plotOutput("fit", height = "430px"))
),
## details of the iterations, shown under the plots
hr(),
fluidRow(
column(7, tableOutput("tab"), htmlOutput("note")),
column(5, helpText("Try the starting values (0, 3) and then (0, 5), first ",
"with both safeguards off, then with each in turn. From ",
"(0, 5) the plain iteration stops because the ",
"information matrix goes numerically singular; the ",
"robustification is what keeps it invertible."))
)
)
### ---------- server ---------------------------------------------------------
server <- function(input, output, session) {
playing <- reactiveVal(FALSE)
play_ms <- 900
## clicking the contour plot moves the starting value
observeEvent(input$click, {
p <- c(input$click$x, input$click$y)
if (all(is.finite(p))) {
updateNumericInput(session, "b0", value = round(p[1], 2))
updateNumericInput(session, "b1", value = round(p[2], 2))
}
})
start <- reactive({
req(is.finite(input$b0), is.finite(input$b1))
c(input$b0, input$b1)
})
modified <- reactive(input$ls != "none" || input$robust != "none")
path <- reactive(nr_path(start(), maxit = 15, ls = input$ls,
robust = input$robust))
## the plain (unsafeguarded) path, shown for comparison
ref <- reactive({
if (!isTRUE(input$compare) || !modified()) return(NULL)
nr_path(start(), maxit = 15)
})
## current iteration, clamped to the length of the path
k <- reactive(min(as.integer(input$iter), path()$n - 1))
## a new path restarts the display at iteration 0
observeEvent(path(), {
updateSliderInput(session, "iter", max = max(1, path()$n - 1), value = 0)
playing(isTRUE(input$autoplay) && path()$n > 1)
})
go_to <- function(i) {
updateSliderInput(session, "iter",
value = max(0, min(i, path()$n - 1)))
}
observeEvent(input$play, {
if (k() >= path()$n - 1) go_to(0)
if (path()$n > 1) playing(TRUE)
})
observeEvent(input$pause, playing(FALSE))
observeEvent(input$step, { playing(FALSE); go_to(k() + 1) })
observeEvent(input$back, { playing(FALSE); go_to(k() - 1) })
observeEvent(input$first, { playing(FALSE); go_to(0) })
observeEvent(input$toend, { playing(FALSE); go_to(path()$n - 1) })
observe({
if (!isTRUE(playing())) return()
kk <- isolate(k())
if (kk >= isolate(path()$n) - 1) { playing(FALSE); return() }
invalidateLater(play_ms, session)
isolate(go_to(kk + 1))
})
output$contour <- renderPlot({
P <- path()$P; kk <- k() + 1L
RP <- ref()
par(mar = c(4, 4, 2, 1))
contour(B0, B1, Z, nlevels = 40, drawlabels = FALSE,
col = "grey70",
xlab = expression(beta[0]), ylab = expression(beta[1]),
main = "negative log-likelihood surface and NR path")
points(0, 1.5, pch = 3, col = "darkgreen", lwd = 2)
## the unsafeguarded path, drawn underneath for comparison
if (!is.null(RP)) {
Q <- RP$P[seq_len(min(kk, RP$n)), , drop = FALSE]
if (nrow(Q) > 1) lines(Q[, 1], Q[, 2], col = "grey45", lwd = 2,
lty = 2)
points(Q[, 1], Q[, 2], pch = 21, bg = "grey85", col = "grey45",
cex = 0.9)
}
if (kk > 1) lines(P[1:kk, 1], P[1:kk, 2], col = "steelblue", lwd = 2)
points(P[1:kk, 1], P[1:kk, 2], pch = 21, bg = "steelblue", cex = 1.1)
cur <- P[kk, ]
inside <- cur[1] >= min(B0) && cur[1] <= max(B0) &&
cur[2] >= min(B1) && cur[2] <= max(B1)
if (inside) {
## Arrow along the score (uphill in the log-likelihood, i.e. minus
## the gradient of negll), of fixed length on screen. The gradient
## is a covector: because the two axes are drawn at different
## scales, its on-screen components are obtained by MULTIPLYING by
## the units-per-inch factors (not dividing, as for a displacement).
## Only then is the arrow perpendicular to the contours as drawn.
g <- -cur[4:5]
upi <- c(diff(par("usr")[1:2]) / par("pin")[1],
diff(par("usr")[3:4]) / par("pin")[2]) # units per inch
gs <- g * upi # inches
## at a stationary point the direction is numerical noise, so the
## arrow is dropped rather than drawn at full length
if (sqrt(sum(g^2)) > 1e-6 && sqrt(sum(gs^2)) > 0) {
d <- 0.45 * gs / sqrt(sum(gs^2)) * upi # 0.45 inch long
arrows(cur[1], cur[2], cur[1] + d[1], cur[2] + d[2],
col = "#E8820C", lwd = 3, length = 0.09)
}
points(cur[1], cur[2], pch = 21, bg = "firebrick", cex = 1.8)
} else {
mtext(sprintf("current iterate (%.3g, %.3g) is off the plotted region",
cur[1], cur[2]), side = 3, line = -1.5,
col = "firebrick", cex = 0.9)
}
leg <- list(txt = c("current iterate", "true (0, 1.5)",
"gradient of log-likelihood (direction)"),
pch = c(21, 3, NA), bg = c("firebrick", NA, NA),
lwd = c(NA, 2, 3), lty = c(NA, NA, 1),
col = c("black", "darkgreen", "#E8820C"))
if (!is.null(RP)) {
leg$txt <- c(leg$txt, "plain Newton-Raphson")
leg$pch <- c(leg$pch, NA); leg$bg <- c(leg$bg, NA)
leg$lwd <- c(leg$lwd, 2); leg$lty <- c(leg$lty, 2)
leg$col <- c(leg$col, "grey45")
}
legend("bottomright", bg = "#ffffffcc", box.col = NA, cex = 0.85,
pch = leg$pch, pt.bg = leg$bg, lwd = leg$lwd, lty = leg$lty,
col = leg$col, legend = leg$txt)
})
## the logistic curve fitted at the current iterate, over the data; the 0/1
## responses are jittered vertically so that their density can be seen
output$fit <- renderPlot({
P <- path()$P; kk <- k() + 1L; RP <- ref()
xs <- seq(min(z), max(z), length.out = 200)
curve_at <- function(b) plogis(b[1] + b[2] * xs)
par(mar = c(4, 4, 2, 1))
plot(z, r_jit, pch = 19, cex = 0.7, col = adjustcolor("grey30", 0.4),
ylim = c(-0.2, 1.2), xlab = "x", ylab = "P(y = 1)",
main = sprintf("fitted logistic curve, iteration %d", kk - 1))
abline(h = c(0, 1), col = "grey85")
lines(xs, curve_at(c(0, 1.5)), col = "darkgreen", lwd = 2, lty = 3)
## the earlier iterates, faint
if (kk > 1) for (i in seq_len(kk - 1))
lines(xs, curve_at(P[i, 1:2]), col = adjustcolor("steelblue", 0.35))
## the plain (unsafeguarded) iterate at the same step
if (!is.null(RP))
lines(xs, curve_at(RP$P[min(kk, RP$n), 1:2]), col = "grey45",
lwd = 2, lty = 2)
lines(xs, curve_at(P[kk, 1:2]), col = "firebrick", lwd = 3)
leg <- list(txt = c("data (jittered)", "current fit", "earlier iterates",
"true curve"),
pch = c(19, NA, NA, NA), lwd = c(NA, 3, 1, 2),
lty = c(NA, 1, 1, 3),
col = c(adjustcolor("grey30", 0.6), "firebrick", "steelblue",
"darkgreen"))
if (!is.null(RP)) {
leg$txt <- c(leg$txt, "plain Newton-Raphson")
leg$pch <- c(leg$pch, NA); leg$lwd <- c(leg$lwd, 2)
leg$lty <- c(leg$lty, 2); leg$col <- c(leg$col, "grey45")
}
legend("right", bg = "#ffffffcc", box.col = NA, cex = 0.85,
pch = leg$pch, lwd = leg$lwd, lty = leg$lty, col = leg$col,
legend = leg$txt)
})
## one decimal throughout; a diverging run reaches values that would
## otherwise be far too wide, so those switch to one-decimal scientific
f1 <- function(x) ifelse(!is.finite(x), "",
ifelse(abs(x) >= 1e5 | (x != 0 & abs(x) < 1e-3),
formatC(x, format = "e", digits = 1),
formatC(x, format = "f", digits = 1)))
## why the run ended, shown as one line under the table
output$note <- renderUI({
R <- path(); RP <- ref()
if (k() + 1 < R$n) return(NULL)
HTML(paste0("<i>", R$note,
if (!is.null(RP)) paste0("; plain Newton-Raphson: ",
RP$note),
"</i>"))
})
output$tab <- renderTable({
P <- path()$P[seq_len(k() + 1), , drop = FALSE]
d <- data.frame(iteration = as.character(seq_len(nrow(P)) - 1),
beta0 = f1(P[, 1]), beta1 = f1(P[, 2]),
negll = f1(P[, 3]), step = f1(P[, 6]),
`cond(J)` = f1(P[, 7]),
check.names = FALSE, stringsAsFactors = FALSE)
tail(d, 5)
}, rownames = FALSE, align = "r", width = "100%")
}
shinyApp(ui, server)This app accompanies Optimization for Maximum Likelihood Estimation in the book.