Regression and Other Stories: Arsenic

Building a logistic regression model: wells in Bangladesh. See Chapter 13 in Regression and Other Stories.

This version uses algorithm='optimizing', which finds the posterior mode, computes Hessian and uses a normal distribution approximation of the posterior. This can work well for non-hierarchical generalized linear models when the number of observations is much higher than the number of covariates.

---

Load packages

library("rprojroot")
root<-has_file(".ROS-Examples-root")$make_fix_file()
library("rstanarm")
library("loo")
invlogit <- plogis

Load data

wells <- read.csv(root("Arsenic/data","wells.csv"))
wells$y <- wells$switch
n <- nrow(wells)

Null model

Log-score for coin flipping

prob <- 0.5
round(log(prob)*sum(wells$y) + log(1-prob)*sum(1-wells$y),1)

[1] -2093.3

Log-score for intercept model

round(prob <- mean(wells$y),2)

[1] 0.58

round(log(prob)*sum(wells$y) + log(1-prob)*sum(1-wells$y),1)

[1] -2059

A single predictor

Fit a model using distance to the nearest safe well

fit_1 <- stan_glm(y ~ dist, family = binomial(link = "logit"), data=wells, algorithm='optimizing')

print(fit_1, digits=3)

stan_glm
 family:       binomial [logit]
 formula:      y ~ dist
 observations: 3020
 predictors:   2
------
            Median MAD_SD
(Intercept)  0.609  0.063
dist        -0.006  0.001

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg

summary(fit_1, digits=3)

Model Info:
 function:     stan_glm
 family:       binomial [logit]
 formula:      y ~ dist
 algorithm:    optimizing
 priors:       see help('prior_summary')
 observations: 3020
 predictors:   2

Estimates:
              Median   MAD_SD   10%    50%    90% 
(Intercept)  0.609    0.063    0.529  0.609  0.689
dist        -0.006    0.001   -0.007 -0.006 -0.005

Monte Carlo diagnostics
            mcse  khat  n_eff
(Intercept) 0.002 0.224 738  
dist        0.000 0.360 752  

For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and khat is the Pareto k diagnostic for importance sampling (perfomance is usually good when khat < 0.7).

LOO log score

(loo1 <- loo(fit_1))

Computed from 1000 by 3020 log-likelihood matrix

         Estimate   SE
elpd_loo  -2040.1 10.4
p_loo         2.0  0.0
looic      4080.3 20.9
------
Monte Carlo SE of elpd_loo is 0.1.

All Pareto k estimates are good (k < 0.5).
See help('pareto-k-diagnostic') for details.

Histogram of distances

hist(wells$dist, breaks=seq(0,10+max(wells$dist),10), freq=TRUE,
     xlab="Distance (in meters) to nearest safe well", ylab="", main="", mgp=c(2,.5,0))

Scale distance in meters to distance in 100 meters

wells$dist100 <- wells$dist/100

Fit a model using scaled distance to the nearest safe well

fit_2 <- stan_glm(y ~ dist100, family = binomial(link = "logit"), data=wells, algorithm='optimizing')

print(fit_2, digits=2)

stan_glm
 family:       binomial [logit]
 formula:      y ~ dist100
 observations: 3020
 predictors:   2
------
            Median MAD_SD
(Intercept)  0.61   0.07 
dist100     -0.63   0.11 

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg

summary(fit_2, digits=2)

Model Info:
 function:     stan_glm
 family:       binomial [logit]
 formula:      y ~ dist100
 algorithm:    optimizing
 priors:       see help('prior_summary')
 observations: 3020
 predictors:   2

Estimates:
              Median   MAD_SD   10%   50%   90%
(Intercept)  0.61     0.07     0.53  0.61  0.69
dist100     -0.63     0.11    -0.76 -0.63 -0.48

Monte Carlo diagnostics
            mcse khat n_eff
(Intercept) 0.00 0.18 761  
dist100     0.00 0.17 717  

For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and khat is the Pareto k diagnostic for importance sampling (perfomance is usually good when khat < 0.7).

LOO log score

(loo2 <- loo(fit_2, save_psis = TRUE))

Computed from 1000 by 3020 log-likelihood matrix

         Estimate   SE
elpd_loo  -2040.3 10.4
p_loo         2.1  0.1
looic      4080.5 20.8
------
Monte Carlo SE of elpd_loo is 0.1.

All Pareto k estimates are good (k < 0.5).
See help('pareto-k-diagnostic') for details.

Plot model fit

jitter_binary <- function(a, jitt=.05){
  a + (1-2*a)*runif(length(a),0,jitt)
}

plot(c(0,max(wells$dist, na.rm=TRUE)*1.02), c(0,1),
     xlab="Distance (in meters) to nearest safe well", ylab="Pr (switching)",
     type="n", xaxs="i", yaxs="i", mgp=c(2,.5,0))
curve(invlogit(coef(fit_1)[1]+coef(fit_1)[2]*x), lwd=1, add=TRUE)
points(wells$dist, jitter_binary(wells$y), pch=20, cex=.1)

Plot uncertainty in the estimated coefficients

sims <- as.matrix(fit_2)
par(pty="s")
plot(sims[1:500,1], sims[1:500,2], xlim=c(.4,.8), ylim=c(-1,0),
     xlab=expression(beta[0]), ylab=expression(beta[1]), mgp=c(1.5,.5,0),
     pch=20, cex=.5, xaxt="n", yaxt="n")
axis(1, seq(.4,.8,.2), mgp=c(1.5,.5,0))
axis(2, seq(-1,0,.5), mgp=c(1.5,.5,0))

Plot uncertainty in the estimated predictions

plot(c(0,max(wells$dist, na.rm=T)*1.02), c(0,1),
     xlab="Distance (in meters) to nearest safe well", ylab="Pr (switching)",
     type="n", xaxs="i", yaxs="i", mgp=c(2,.5,0))
for (j in 1:20) {
    curve (invlogit(sims[j,1]+sims[j,2]*x/100), lwd=.5,
           col="darkgray", add=TRUE)
}
curve(invlogit(coef(fit_2)[1]+coef(fit_2)[2]*x/100), lwd=1, add=T)
points(wells$dist, jitter_binary(wells$y), pch=20, cex=.1)

Two predictors

Histogram of arsenic levels

hist(wells$arsenic, breaks=seq(0,.25+max(wells$arsenic),.25), freq=TRUE,
     xlab="Arsenic concentration in well water", ylab="", main="", mgp=c(2,.5,0))

Fit a model using scaled distance and arsenic level

fit_3 <- stan_glm(y ~ dist100 + arsenic, family = binomial(link = "logit"),
                  data=wells, algorithm='optimizing')

print(fit_3, digits=2)

stan_glm
 family:       binomial [logit]
 formula:      y ~ dist100 + arsenic
 observations: 3020
 predictors:   3
------
            Median MAD_SD
(Intercept)  0.01   0.08 
dist100     -0.90   0.11 
arsenic      0.46   0.04 

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg

summary(fit_3, digits=2)

Model Info:
 function:     stan_glm
 family:       binomial [logit]
 formula:      y ~ dist100 + arsenic
 algorithm:    optimizing
 priors:       see help('prior_summary')
 observations: 3020
 predictors:   3

Estimates:
              Median   MAD_SD   10%   50%   90%
(Intercept)  0.01     0.08    -0.10  0.01  0.10
dist100     -0.90     0.11    -1.03 -0.90 -0.76
arsenic      0.46     0.04     0.41  0.46  0.51

Monte Carlo diagnostics
            mcse khat n_eff
(Intercept) 0.00 0.24 704  
dist100     0.00 0.14 762  
arsenic     0.00 0.20 765  

For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and khat is the Pareto k diagnostic for importance sampling (perfomance is usually good when khat < 0.7).

LOO log score

(loo3 <- loo(fit_3, save_psis = TRUE))

Computed from 1000 by 3020 log-likelihood matrix

         Estimate   SE
elpd_loo  -1968.4 15.6
p_loo         3.1  0.1
looic      3936.7 31.3
------
Monte Carlo SE of elpd_loo is 0.1.

All Pareto k estimates are good (k < 0.5).
See help('pareto-k-diagnostic') for details.

Compare models

loo_compare(loo2, loo3)

      elpd_diff se_diff
fit_3   0.0       0.0  
fit_2 -71.9      12.1  

Average improvement in LOO predictive probabilities

from dist100 to dist100 + arsenic

pred2 <- loo_predict(fit_2, psis_object = loo2$psis_object)$value
pred3 <- loo_predict(fit_3, psis_object = loo3$psis_object)$value
round(mean(c(pred3[wells$y==1]-pred2[wells$y==1],pred2[wells$y==0]-pred3[wells$y==0])),3)

[1] 0.022

Plot model fits

plot(c(0,max(wells$dist,na.rm=T)*1.02), c(0,1),
     xlab="Distance (in meters) to nearest safe well", ylab="Pr (switching)",
     type="n", xaxs="i", yaxs="i", mgp=c(2,.5,0))
points(wells$dist, jitter_binary(wells$y), pch=20, cex=.1)
curve(invlogit(coef(fit_3)[1]+coef(fit_3)[2]*x/100+coef(fit_3)[3]*.50), lwd=.5, add=T)
curve(invlogit(coef(fit_3)[1]+coef(fit_3)[2]*x/100+coef(fit_3)[3]*1.00), lwd=.5, add=T)
text(50, .27, "if As = 0.5", adj=0, cex=.8)
text(75, .50, "if As = 1.0", adj=0, cex=.8)

plot(c(0,max(wells$arsenic,na.rm=T)*1.02), c(0,1),
     xlab="Arsenic concentration in well water", ylab="Pr (switching)",
     type="n", xaxs="i", yaxs="i", mgp=c(2,.5,0))
points(wells$arsenic, jitter_binary(wells$y), pch=20, cex=.1)
curve(invlogit(coef(fit_3)[1]+coef(fit_3)[2]*0+coef(fit_3)[3]*x), from=0.5, lwd=.5, add=T)
curve(invlogit(coef(fit_3)[1]+coef(fit_3)[2]*0.5+coef(fit_3)[3]*x), from=0.5, lwd=.5, add=T)
text(.5, .78, "if dist = 0", adj=0, cex=.8)
text(2, .6, "if dist = 50", adj=0, cex=.8)

Interaction

Fit a model using scaled distance, arsenic level, and an interaction

fit_4 <- stan_glm(y ~ dist100 + arsenic + dist100:arsenic,
                  family = binomial(link="logit"), data = wells, algorithm='optimizing')

print(fit_4, digits=2)

stan_glm
 family:       binomial [logit]
 formula:      y ~ dist100 + arsenic + dist100:arsenic
 observations: 3020
 predictors:   4
------
                Median MAD_SD
(Intercept)     -0.14   0.12 
dist100         -0.59   0.21 
arsenic          0.56   0.07 
dist100:arsenic -0.17   0.10 

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg

summary(fit_4, digits=2)

Model Info:
 function:     stan_glm
 family:       binomial [logit]
 formula:      y ~ dist100 + arsenic + dist100:arsenic
 algorithm:    optimizing
 priors:       see help('prior_summary')
 observations: 3020
 predictors:   4

Estimates:
                  Median   MAD_SD   10%   50%   90%
(Intercept)     -0.14     0.12    -0.28 -0.14  0.01
dist100         -0.59     0.21    -0.86 -0.59 -0.31
arsenic          0.56     0.07     0.47  0.56  0.65
dist100:arsenic -0.17     0.10    -0.31 -0.17 -0.06

Monte Carlo diagnostics
                mcse khat n_eff
(Intercept)     0.00 0.22 762  
dist100         0.01 0.39 745  
arsenic         0.00 0.20 778  
dist100:arsenic 0.00 0.28 762  

For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and khat is the Pareto k diagnostic for importance sampling (perfomance is usually good when khat < 0.7).

LOO log score

(loo4 <- loo(fit_4))

Computed from 1000 by 3020 log-likelihood matrix

         Estimate   SE
elpd_loo  -1968.1 15.9
p_loo         4.5  0.3
looic      3936.2 31.8
------
Monte Carlo SE of elpd_loo is 0.1.

All Pareto k estimates are good (k < 0.5).
See help('pareto-k-diagnostic') for details.

Compare models

loo_compare(loo3, loo4)

      elpd_diff se_diff
fit_4  0.0       0.0   
fit_3 -0.2       1.9   

Centering the input variables

wells$c_dist100 <- wells$dist100 - mean(wells$dist100)
wells$c_arsenic <- wells$arsenic - mean(wells$arsenic)
fit_5 <- stan_glm(y ~ c_dist100 + c_arsenic + c_dist100:c_arsenic,
                  family = binomial(link="logit"), data = wells, algorithm='optimizing')

print(fit_5, digits=2)

stan_glm
 family:       binomial [logit]
 formula:      y ~ c_dist100 + c_arsenic + c_dist100:c_arsenic
 observations: 3020
 predictors:   4
------
                    Median MAD_SD
(Intercept)          0.35   0.04 
c_dist100           -0.88   0.10 
c_arsenic            0.47   0.04 
c_dist100:c_arsenic -0.17   0.10 

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg

summary(fit_5, digits=2)

Model Info:
 function:     stan_glm
 family:       binomial [logit]
 formula:      y ~ c_dist100 + c_arsenic + c_dist100:c_arsenic
 algorithm:    optimizing
 priors:       see help('prior_summary')
 observations: 3020
 predictors:   4

Estimates:
                      Median   MAD_SD   10%   50%   90%
(Intercept)          0.35     0.04     0.30  0.35  0.40
c_dist100           -0.88     0.10    -1.01 -0.88 -0.75
c_arsenic            0.47     0.04     0.42  0.47  0.52
c_dist100:c_arsenic -0.17     0.10    -0.31 -0.17 -0.04

Monte Carlo diagnostics
                    mcse khat n_eff
(Intercept)         0.00 0.26 776  
c_dist100           0.00 0.19 818  
c_arsenic           0.00 0.21 785  
c_dist100:c_arsenic 0.00 0.20 779  

For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and khat is the Pareto k diagnostic for importance sampling (perfomance is usually good when khat < 0.7).

Plot model fits

plot(c(0,max(wells$dist,na.rm=T)*1.02), c(0,1),
     xlab="Distance (in meters) to nearest safe well", ylab="Pr (switching)",
     type="n", xaxs="i", yaxs="i", mgp=c(2,.5,0))
points(wells$dist, jitter_binary(wells$y), pch=20, cex=.1)
curve(invlogit(coef(fit_4)[1]+coef(fit_4)[2]*x/100+coef(fit_4)[3]*.50+coef(fit_4)[4]*x/100*.50), lwd=.5, add=T)
curve(invlogit(coef(fit_4)[1]+coef(fit_4)[2]*x/100+coef(fit_4)[3]*1.00+coef(fit_4)[4]*x/100*1.00), lwd=.5, add=T)
text (50, .29, "if As = 0.5", adj=0, cex=.8)
text (75, .50, "if As = 1.0", adj=0, cex=.8)

plot(c(0,max(wells$arsenic,na.rm=T)*1.02), c(0,1),
     xlab="Arsenic concentration in well water", ylab="Pr (switching)",
     type="n", xaxs="i", yaxs="i", mgp=c(2,.5,0))
points(wells$arsenic, jitter_binary(wells$y), pch=20, cex=.1)
curve(invlogit(coef(fit_4)[1]+coef(fit_4)[2]*0+coef(fit_4)[3]*x+coef(fit_4)[4]*0*x), from=0.5, lwd=.5, add=T)
curve(invlogit(coef(fit_4)[1]+coef(fit_4)[2]*0.5+coef(fit_4)[3]*x+coef(fit_4)[4]*0.5*x), from=0.5, lwd=.5, add=T)
text (.5, .78, "if dist = 0", adj=0, cex=.8)
text (2, .6, "if dist = 50", adj=0, cex=.8)

More predictors

Adding social predictors

fit_6 <- stan_glm(y ~ dist100 + arsenic + educ4 + assoc,
                  family = binomial(link="logit"), data = wells, algorithm='optimizing')

print(fit_6, digits=2)

stan_glm
 family:       binomial [logit]
 formula:      y ~ dist100 + arsenic + educ4 + assoc
 observations: 3020
 predictors:   5
------
            Median MAD_SD
(Intercept) -0.16   0.11 
dist100     -0.90   0.10 
arsenic      0.47   0.04 
educ4        0.17   0.04 
assoc       -0.13   0.07 

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg

summary(fit_6, digits=2)

Model Info:
 function:     stan_glm
 family:       binomial [logit]
 formula:      y ~ dist100 + arsenic + educ4 + assoc
 algorithm:    optimizing
 priors:       see help('prior_summary')
 observations: 3020
 predictors:   5

Estimates:
              Median   MAD_SD   10%   50%   90%
(Intercept) -0.16     0.11    -0.29 -0.16 -0.03
dist100     -0.90     0.10    -1.03 -0.90 -0.77
arsenic      0.47     0.04     0.42  0.47  0.52
educ4        0.17     0.04     0.12  0.17  0.22
assoc       -0.13     0.07    -0.22 -0.13 -0.03

Monte Carlo diagnostics
            mcse khat n_eff
(Intercept) 0.00 0.26 719  
dist100     0.00 0.28 732  
arsenic     0.00 0.27 738  
educ4       0.00 0.25 734  
assoc       0.00 0.26 722  

For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and khat is the Pareto k diagnostic for importance sampling (perfomance is usually good when khat < 0.7).

LOO log score

(loo6 <- loo(fit_6))

Computed from 1000 by 3020 log-likelihood matrix

         Estimate   SE
elpd_loo  -1959.0 16.1
p_loo         5.2  0.1
looic      3918.0 32.3
------
Monte Carlo SE of elpd_loo is 0.1.

All Pareto k estimates are good (k < 0.5).
See help('pareto-k-diagnostic') for details.

Compare models

loo_compare(loo4, loo6)

      elpd_diff se_diff
fit_6  0.0       0.0   
fit_4 -9.1       5.1   

Remove assoc

fit_7 <- stan_glm(y ~ dist100 + arsenic + educ4,
                  family = binomial(link="logit"), data = wells, algorithm='optimizing')

print(fit_7, digits=2)

stan_glm
 family:       binomial [logit]
 formula:      y ~ dist100 + arsenic + educ4
 observations: 3020
 predictors:   4
------
            Median MAD_SD
(Intercept) -0.21   0.09 
dist100     -0.90   0.10 
arsenic      0.47   0.04 
educ4        0.17   0.04 

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg

summary(fit_7, digits=2)

Model Info:
 function:     stan_glm
 family:       binomial [logit]
 formula:      y ~ dist100 + arsenic + educ4
 algorithm:    optimizing
 priors:       see help('prior_summary')
 observations: 3020
 predictors:   4

Estimates:
              Median   MAD_SD   10%   50%   90%
(Intercept) -0.21     0.09    -0.33 -0.21 -0.09
dist100     -0.90     0.10    -1.02 -0.90 -0.76
arsenic      0.47     0.04     0.42  0.47  0.52
educ4        0.17     0.04     0.12  0.17  0.22

Monte Carlo diagnostics
            mcse khat n_eff
(Intercept) 0.00 0.31 775  
dist100     0.00 0.27 768  
arsenic     0.00 0.40 790  
educ4       0.00 0.36 790  

For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and khat is the Pareto k diagnostic for importance sampling (perfomance is usually good when khat < 0.7).

LOO log score

(loo7 <- loo(fit_7))

Computed from 1000 by 3020 log-likelihood matrix

         Estimate   SE
elpd_loo  -1959.4 16.0
p_loo         4.2  0.1
looic      3918.7 32.0
------
Monte Carlo SE of elpd_loo is 0.1.

All Pareto k estimates are good (k < 0.5).
See help('pareto-k-diagnostic') for details.

Compare models

loo_compare(loo4, loo7)

      elpd_diff se_diff
fit_7  0.0       0.0   
fit_4 -8.7       4.8   

loo_compare(loo6, loo7)

      elpd_diff se_diff
fit_6  0.0       0.0   
fit_7 -0.4       1.6   

Add interactions with education

wells$c_educ4 <- wells$educ4 - mean(wells$educ4)
fit_8 <- stan_glm(y ~ c_dist100 + c_arsenic + c_educ4 +
                      c_dist100:c_educ4 + c_arsenic:c_educ4,
                  family = binomial(link="logit"), data = wells, algorithm='optimizing')

print(fit_8, digits=2)

stan_glm
 family:       binomial [logit]
 formula:      y ~ c_dist100 + c_arsenic + c_educ4 + c_dist100:c_educ4 + c_arsenic:c_educ4
 observations: 3020
 predictors:   6
------
                  Median MAD_SD
(Intercept)        0.35   0.04 
c_dist100         -0.92   0.11 
c_arsenic          0.49   0.04 
c_educ4            0.19   0.04 
c_dist100:c_educ4  0.33   0.11 
c_arsenic:c_educ4  0.08   0.04 

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg

summary(fit_8, digits=2)

Model Info:
 function:     stan_glm
 family:       binomial [logit]
 formula:      y ~ c_dist100 + c_arsenic + c_educ4 + c_dist100:c_educ4 + c_arsenic:c_educ4
 algorithm:    optimizing
 priors:       see help('prior_summary')
 observations: 3020
 predictors:   6

Estimates:
                    Median   MAD_SD   10%   50%   90%
(Intercept)        0.35     0.04     0.29  0.35  0.40
c_dist100         -0.92     0.11    -1.06 -0.92 -0.78
c_arsenic          0.49     0.04     0.44  0.49  0.54
c_educ4            0.19     0.04     0.14  0.19  0.24
c_dist100:c_educ4  0.33     0.11     0.20  0.33  0.48
c_arsenic:c_educ4  0.08     0.04     0.02  0.08  0.13

Monte Carlo diagnostics
                  mcse khat n_eff
(Intercept)       0.00 0.37 726  
c_dist100         0.00 0.54 689  
c_arsenic         0.00 0.40 743  
c_educ4           0.00 0.52 780  
c_dist100:c_educ4 0.00 0.36 714  
c_arsenic:c_educ4 0.00 0.44 732  

For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and khat is the Pareto k diagnostic for importance sampling (perfomance is usually good when khat < 0.7).

LOO log score

(loo8 <- loo(fit_8, save_psis=TRUE))

Computed from 1000 by 3020 log-likelihood matrix

         Estimate   SE
elpd_loo  -1952.8 16.5
p_loo         6.5  0.3
looic      3905.7 33.0
------
Monte Carlo SE of elpd_loo is 0.1.

All Pareto k estimates are good (k < 0.5).
See help('pareto-k-diagnostic') for details.

Compare models

loo_compare(loo3, loo8)

      elpd_diff se_diff
fit_8   0.0       0.0  
fit_3 -15.5       6.3  

loo_compare(loo7, loo8)

      elpd_diff se_diff
fit_8  0.0       0.0   
fit_7 -6.5       4.3   

Average improvement in LOO predictive probabilities

from dist100 + arsenic to dist100 + arsenic + educ4 + dist100:educ4 + arsenic:educ4

pred8 <- loo_predict(fit_8, psis_object = loo8$psis_object)$value
round(mean(c(pred8[wells$y==1]-pred3[wells$y==1],pred3[wells$y==0]-pred8[wells$y==0])),3)

[1] 0.005

Transformation of variable

Fit a model using scaled distance and log arsenic level

wells$log_arsenic <- log(wells$arsenic)

fit_3a <- stan_glm(y ~ dist100 + log_arsenic, family = binomial(link = "logit"),
                   data = wells, algorithm='optimizing')

print(fit_3a, digits=2)

stan_glm
 family:       binomial [logit]
 formula:      y ~ dist100 + log_arsenic
 observations: 3020
 predictors:   3
------
            Median MAD_SD
(Intercept)  0.52   0.06 
dist100     -0.97   0.10 
log_arsenic  0.88   0.07 

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg

summary(fit_3a, digits=2)

Model Info:
 function:     stan_glm
 family:       binomial [logit]
 formula:      y ~ dist100 + log_arsenic
 algorithm:    optimizing
 priors:       see help('prior_summary')
 observations: 3020
 predictors:   3

Estimates:
              Median   MAD_SD   10%   50%   90%
(Intercept)  0.52     0.06     0.45  0.52  0.61
dist100     -0.97     0.10    -1.12 -0.97 -0.85
log_arsenic  0.88     0.07     0.79  0.88  0.97

Monte Carlo diagnostics
            mcse khat n_eff
(Intercept) 0.00 0.02 782  
dist100     0.00 0.07 801  
log_arsenic 0.00 0.01 779  

For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and khat is the Pareto k diagnostic for importance sampling (perfomance is usually good when khat < 0.7).

LOO log score

(loo3a <- loo(fit_3a))

Computed from 1000 by 3020 log-likelihood matrix

         Estimate   SE
elpd_loo  -1952.3 16.3
p_loo         3.1  0.1
looic      3904.6 32.7
------
Monte Carlo SE of elpd_loo is 0.1.

All Pareto k estimates are good (k < 0.5).
See help('pareto-k-diagnostic') for details.

Compare models

loo_compare(loo3, loo3a)

       elpd_diff se_diff
fit_3a   0.0       0.0  
fit_3  -16.1       4.4  

Fit a model using scaled distance, log arsenic level, and an interaction

fit_4a <- stan_glm(y ~ dist100 + log_arsenic + dist100:log_arsenic,
                  family = binomial(link = "logit"), data = wells, algorithm='optimizing')

print(fit_4a, digits=2)

stan_glm
 family:       binomial [logit]
 formula:      y ~ dist100 + log_arsenic + dist100:log_arsenic
 observations: 3020
 predictors:   4
------
                    Median MAD_SD
(Intercept)          0.49   0.07 
dist100             -0.87   0.13 
log_arsenic          0.98   0.10 
dist100:log_arsenic -0.23   0.17 

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg

summary(fit_4a, digits=2)

Model Info:
 function:     stan_glm
 family:       binomial [logit]
 formula:      y ~ dist100 + log_arsenic + dist100:log_arsenic
 algorithm:    optimizing
 priors:       see help('prior_summary')
 observations: 3020
 predictors:   4

Estimates:
                      Median   MAD_SD   10%   50%   90%
(Intercept)          0.49     0.07     0.41  0.49  0.58
dist100             -0.87     0.13    -1.06 -0.87 -0.70
log_arsenic          0.98     0.10     0.85  0.98  1.12
dist100:log_arsenic -0.23     0.17    -0.46 -0.23  0.00

Monte Carlo diagnostics
                    mcse khat n_eff
(Intercept)         0.00 0.29 758  
dist100             0.00 0.27 767  
log_arsenic         0.00 0.21 714  
dist100:log_arsenic 0.01 0.28 750  

For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and khat is the Pareto k diagnostic for importance sampling (perfomance is usually good when khat < 0.7).

LOO log score

(loo4a <- loo(fit_4a))

Computed from 1000 by 3020 log-likelihood matrix

         Estimate   SE
elpd_loo  -1952.5 16.4
p_loo         4.1  0.1
looic      3904.9 32.8
------
Monte Carlo SE of elpd_loo is 0.1.

All Pareto k estimates are good (k < 0.5).
See help('pareto-k-diagnostic') for details.

Compare models

loo_compare(loo3a, loo4a)

       elpd_diff se_diff
fit_3a  0.0       0.0   
fit_4a -0.2       1.2   

Add interactions with education

wells$c_log_arsenic <- wells$log_arsenic - mean(wells$log_arsenic)

fit_8a <- stan_glm(y ~ c_dist100 + c_log_arsenic + c_educ4 +
                      c_dist100:c_educ4 + c_log_arsenic:c_educ4,
                  family = binomial(link="logit"), data = wells, algorithm='optimizing')

print(fit_8a, digits=2)

stan_glm
 family:       binomial [logit]
 formula:      y ~ c_dist100 + c_log_arsenic + c_educ4 + c_dist100:c_educ4 + 
       c_log_arsenic:c_educ4
 observations: 3020
 predictors:   6
------
                      Median MAD_SD
(Intercept)            0.34   0.04 
c_dist100             -1.01   0.11 
c_log_arsenic          0.91   0.07 
c_educ4                0.18   0.04 
c_dist100:c_educ4      0.35   0.11 
c_log_arsenic:c_educ4  0.07   0.07 

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg

summary(fit_8a, digits=2)

Model Info:
 function:     stan_glm
 family:       binomial [logit]
 formula:      y ~ c_dist100 + c_log_arsenic + c_educ4 + c_dist100:c_educ4 + 
       c_log_arsenic:c_educ4
 algorithm:    optimizing
 priors:       see help('prior_summary')
 observations: 3020
 predictors:   6

Estimates:
                        Median   MAD_SD   10%   50%   90%
(Intercept)            0.34     0.04     0.28  0.34  0.38
c_dist100             -1.01     0.11    -1.15 -1.01 -0.88
c_log_arsenic          0.91     0.07     0.82  0.91  1.00
c_educ4                0.18     0.04     0.13  0.18  0.23
c_dist100:c_educ4      0.35     0.11     0.20  0.35  0.49
c_log_arsenic:c_educ4  0.07     0.07    -0.03  0.07  0.15

Monte Carlo diagnostics
                      mcse khat n_eff
(Intercept)           0.00 0.31 700  
c_dist100             0.00 0.29 758  
c_log_arsenic         0.00 0.26 717  
c_educ4               0.00 0.31 736  
c_dist100:c_educ4     0.00 0.30 705  
c_log_arsenic:c_educ4 0.00 0.25 765  

For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and khat is the Pareto k diagnostic for importance sampling (perfomance is usually good when khat < 0.7).

LOO log score

(loo8a <- loo(fit_8a, save_psis=TRUE))

Computed from 1000 by 3020 log-likelihood matrix

         Estimate   SE
elpd_loo  -1938.0 17.2
p_loo         6.2  0.2
looic      3876.1 34.3
------
Monte Carlo SE of elpd_loo is 0.1.

All Pareto k estimates are good (k < 0.5).
See help('pareto-k-diagnostic') for details.

Compare models

loo_compare(loo8, loo8a)

       elpd_diff se_diff
fit_8a   0.0       0.0  
fit_8  -14.8       4.3