Regression and Other Stories: Golf

Gold putting accuracy: Fitting a nonlinear model using Stan. See Chapter 22 in Regression and Other Stories.

---

Load packages

library("rprojroot")
root<-has_file(".ROS-Examples-root")$make_fix_file()
library("rstan")
rstan_options(auto_write = TRUE)
library("rstanarm")
invlogit <- plogis

Set up the data

golf <- read.table(root("Golf/data","golf.txt"), header=TRUE, skip=2)
golf$se <- with(golf, sqrt((y/n)*(1-y/n)/n))
r <- (1.68/2)/12
R <- (4.25/2)/12
golf_data <- list(x=golf$x, y=golf$y, n=golf$n, J=nrow(golf), r=r, R=R)

Plot data

par(mar=c(3,3,2,1), mgp=c(1.7,.5,0), tck=-.02)
with(golf, {
    plot(x, y/n, xlim=c(0, 1.1*max(x)), ylim=c(0, 1.02),
         xaxs="i", yaxs="i", pch=20, bty="l",
         xlab="Distance from hole (feet)", ylab="Probability of success",
         main="Data on putts in pro golf")
    segments(x, y/n + se, x, y/n-se, lwd=.5)
    text(x + .4, y/n + se + .02, paste(y, "/", n,sep=""), cex=.6, col="gray40")
})

Logistic regression model with rstanarm

fit1 <- stan_glm(cbind(y, n-y) ~ x, family=binomial(link="logit"), data=golf,
                 refresh = 0)

Post-processing

a_hat <- fit1$coefficients[1]
b_hat <- fit1$coefficients[2]

Plot logistic regression result

par(mar=c(3,3,2,1), mgp=c(1.7,.5,0), tck=-.02)
with(golf, {
    plot(x, y/n, xlim=c(0, 1.1*max(x)), ylim=c(0, 1.02),
         xaxs="i", yaxs="i", pch=20, bty="l",
         xlab="Distance from hole (feet)", ylab="Probability of success",
         main="Fitted logistic regression")
    segments(x, y/n + se, x, y/n-se, lwd=.5)
    curve(invlogit(a_hat + b_hat*x), from=0, to=1.1*max(x), add=TRUE)
    text(10.6, .57, paste("Logistic regression,\n    a = ",
                          round(a_hat, 2), ", b = ", round(b_hat, 2), sep=""))
})

Logistic regression model with rstan

stanfile_golf_logistic <- root("Golf","golf_logistic.stan")
writeLines(readLines(stanfile_golf_logistic))

data {
  int J;
  int n[J];
  vector[J] x;
  int y[J];
}
parameters {
  real a;
  real b;
}
model {
  y ~ binomial_logit(n, a + b*x);
}

fit_logistic <- stan(file = stanfile_golf_logistic, data = golf_data,
                     refresh = 0)
print(fit_logistic)

Inference for Stan model: golf_logistic.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

         mean se_mean   sd     2.5%      25%      50%      75%    97.5% n_eff
a        2.23    0.00 0.06     2.12     2.19     2.23     2.27     2.35  1340
b       -0.26    0.00 0.01    -0.27    -0.26    -0.26    -0.25    -0.24  1341
lp__ -3021.14    0.03 1.00 -3023.93 -3021.51 -3020.82 -3020.44 -3020.18  1400
     Rhat
a       1
b       1
lp__    1

Samples were drawn using NUTS(diag_e) at Sat Aug 22 12:03:46 2020.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

Post-processing

sims_logistic <- as.matrix(fit_logistic)
a_hat <- median(sims_logistic[,"a"])
b_hat <- median(sims_logistic[,"b"])

Plot logistic regression result

The result is indistinguishable from rstanarm logistic model.

par(mar=c(3,3,2,1), mgp=c(1.7,.5,0), tck=-.02)
with(golf, {
    plot(x, y/n, xlim=c(0, 1.1*max(x)), ylim=c(0, 1.02),
         xaxs="i", yaxs="i", pch=20, bty="l",
         xlab="Distance from hole (feet)", ylab="Probability of success",
         main="Fitted logistic regression")
    segments(x, y/n + se, x, y/n-se, lwd=.5)
    curve(invlogit(a_hat + b_hat*x), from=0, to=1.1*max(x), add=TRUE)
    text(10.6, .57, paste("Logistic regression,\n    a = ",
                          round(a_hat, 2), ", b = ", round(b_hat, 2), sep=""))
})

Geometry based nonlinear model

stanfile_golf_trig <- root("Golf","golf_trig.stan")
writeLines(readLines(stanfile_golf_trig))

data {
  int J;
  int n[J];
  vector[J] x;
  int y[J];
  real r;
  real R;
}
parameters {
  real<lower=0> sigma;
}
model {
  vector[J] p = 2*Phi(asin((R-r) ./ x) / sigma) - 1;
  y ~ binomial(n, p);
}
generated quantities {
  real sigma_degrees;
  sigma_degrees = (180/pi())*sigma;
}

fit_trig <- stan(file = stanfile_golf_trig, data = golf_data, refresh = 0)
print(fit_trig)

Inference for Stan model: golf_trig.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

                  mean se_mean   sd     2.5%      25%      50%      75%
sigma             0.03    0.00 0.00     0.03     0.03     0.03     0.03
sigma_degrees     1.53    0.00 0.02     1.48     1.51     1.53     1.54
lp__          -2926.77    0.02 0.72 -2928.77 -2926.92 -2926.50 -2926.32
                 97.5% n_eff Rhat
sigma             0.03  1522    1
sigma_degrees     1.57  1522    1
lp__          -2926.27  1844    1

Samples were drawn using NUTS(diag_e) at Sat Aug 22 12:04:32 2020.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

Post-processing

sims_trig <- as.matrix(fit_trig)
sigma_hat <- median(sims_trig[,"sigma"])

Plot geometry based model result

par(mar=c(3,3,2,1), mgp=c(1.7,.5,0), tck=-.02)
with(golf, {
    plot(x, y/n, xlim=c(0, 1.1*max(x)), ylim=c(0, 1.02),
         xaxs="i", yaxs="i", pch=20, bty="l",
         xlab="Distance from hole (feet)", ylab="Probability of success",
         main="Custom nonlinear model fit in Stan")
    segments(x, y/n + se, x, y/n-se, lwd=.5)
    x_grid <- seq(R-r, 1.1*max(x), .01)
    p_grid <- 2*pnorm(asin((R-r)/x_grid) / sigma_hat) - 1
    lines(c(0, R-r, x_grid), c(1, 1, p_grid))
    text(18.5, .26, paste("Geometry-based model,\n sigma = ",
                          round(sigma_hat*180/pi, 1), " degrees", sep=""))
})

Plot geometry based model posterior draws of sigma

par(mar=c(3,3,2,1), mgp=c(1.7,.5,0), tck=-.02)
with(golf, {
    plot(x, y/n, xlim=c(0, 1.1*max(x)), ylim=c(0, 1.02),
         xaxs="i", yaxs="i", pch=20, bty="l",
         xlab="Distance from hole (feet)", ylab="Probability of success",
         main="Custom nonlinear model fit in Stan")
    segments(x, y/n + se, x, y/n-se, lwd=.5)
    x_grid <- seq(R-r, 1.1*max(x), .01)
    n_sims <- length(sims_trig[,"sigma"])
    for (s in sample(n_sims, 20)){
        p_grid <- 2*pnorm(asin((R-r)/x_grid) / sims_trig[s,"sigma"]) - 1
        lines(c(0, R-r, x_grid), c(1, 1, p_grid), lwd=0.5)
    }
    text(18.5, .26, "Geometry-based model,\n post draws of sigma")
})

Plot two models in same figure

par(mar=c(3,3,2,1), mgp=c(1.7,.5,0), tck=-.02)
with(golf, {
    plot(x, y/n, xlim=c(0, 1.1*max(x)), ylim=c(0, 1.02),
         xaxs="i", yaxs="i", pch=20, bty="l",
         xlab="Distance from hole (feet)", ylab="Probability of success",
         main="Two models fit to the golf putting data")
    segments(x, y/n + se, x, y/n-se, lwd=.5)
    curve(invlogit(a_hat + b_hat*x), from=0, to=1.1*max(x), add=TRUE)
    x_grid <- seq(R-r, 1.1*max(x), .01)
    p_grid <- 2*pnorm(asin((R-r)/x_grid) / sigma_hat) - 1
    lines(c(0, R-r, x_grid), c(1, 1, p_grid))
    text(10.3, .58, "Logistic regression")
    text(18.5, .24, "Geometry-based model")
})