Regression and Other Stories: Arsenic

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

---

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"))
head(wells)

  switch arsenic   dist dist100 assoc educ educ4
1      1    2.36 16.826 0.16826     0    0  0.00
2      1    0.71 47.322 0.47322     0    0  0.00
3      0    2.07 20.967 0.20967     0   10  2.50
4      1    1.15 21.486 0.21486     0   12  3.00
5      1    1.10 40.874 0.40874     1   14  3.50
6      1    3.90 69.518 0.69518     1    9  2.25

n <- nrow(wells)

Null model

Log-score for coin flipping

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

[1] -2093.3

Log-score for intercept model

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

[1] 0.58

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

[1] -2059

A single predictor

Fit a model using distance to the nearest safe well

fit_1 <- stan_glm(switch ~ dist, family = binomial(link = "logit"), data=wells)

print(fit_1, digits=3)

stan_glm
 family:       binomial [logit]
 formula:      switch ~ dist
 observations: 3020
 predictors:   2
------
            Median MAD_SD
(Intercept)  0.605  0.058
dist        -0.006  0.001

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

LOO log score

(loo1 <- loo(fit_1))

Computed from 4000 by 3020 log-likelihood matrix

         Estimate   SE
elpd_loo  -2040.1 10.4
p_loo         1.9  0.0
looic      4080.1 20.8
------
Monte Carlo SE of elpd_loo is 0.0.

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(switch ~ dist100, family = binomial(link = "logit"), data=wells)

print(fit_2, digits=2)

stan_glm
 family:       binomial [logit]
 formula:      switch ~ dist100
 observations: 3020
 predictors:   2
------
            Median MAD_SD
(Intercept)  0.61   0.06 
dist100     -0.62   0.09 

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

LOO log score

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

Computed from 4000 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.0.

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$switch), 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$switch), 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(switch ~ dist100 + arsenic, family = binomial(link = "logit"),
                  data=wells)

print(fit_3, digits=2)

stan_glm
 family:       binomial [logit]
 formula:      switch ~ 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

LOO log score

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

Computed from 4000 by 3020 log-likelihood matrix

         Estimate   SE
elpd_loo  -1968.5 15.7
p_loo         3.3  0.1
looic      3937.0 31.3
------
Monte Carlo SE of elpd_loo is 0.0.

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.6      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$switch==1]-pred2[wells$switch==1],pred2[wells$switch==0]-pred3[wells$switch==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$switch), 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$switch), 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(switch ~ dist100 + arsenic + dist100:arsenic,
                  family = binomial(link="logit"), data = wells)

print(fit_4, digits=2)

stan_glm
 family:       binomial [logit]
 formula:      switch ~ dist100 + arsenic + dist100:arsenic
 observations: 3020
 predictors:   4
------
                Median MAD_SD
(Intercept)     -0.15   0.11 
dist100         -0.58   0.20 
arsenic          0.55   0.07 
dist100:arsenic -0.18   0.10 

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

LOO log score

(loo4 <- loo(fit_4))

Computed from 4000 by 3020 log-likelihood matrix

         Estimate   SE
elpd_loo  -1968.0 15.9
p_loo         4.3  0.3
looic      3935.9 31.7
------
Monte Carlo SE of elpd_loo is 0.0.

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.5       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(switch ~ c_dist100 + c_arsenic + c_dist100:c_arsenic,
                  family = binomial(link="logit"), data = wells)

print(fit_5, digits=2)

stan_glm
 family:       binomial [logit]
 formula:      switch ~ 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.18   0.10 

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

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$switch), 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$switch), 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(switch ~ dist100 + arsenic + educ4 + assoc,
                  family = binomial(link="logit"), data = wells)

print(fit_6, digits=2)

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

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

LOO log score

(loo6 <- loo(fit_6))

Computed from 4000 by 3020 log-likelihood matrix

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

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.0       5.1   

Remove assoc

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

print(fit_7, digits=2)

stan_glm
 family:       binomial [logit]
 formula:      switch ~ dist100 + arsenic + educ4
 observations: 3020
 predictors:   4
------
            Median MAD_SD
(Intercept) -0.22   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

LOO log score

(loo7 <- loo(fit_7))

Computed from 4000 by 3020 log-likelihood matrix

         Estimate   SE
elpd_loo  -1959.4 16.1
p_loo         4.3  0.1
looic      3918.9 32.1
------
Monte Carlo SE of elpd_loo is 0.0.

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.5       4.8   

loo_compare(loo6, loo7)

      elpd_diff se_diff
fit_6  0.0       0.0   
fit_7 -0.5       1.6   

Add interactions with education

wells$c_educ4 <- wells$educ4 - mean(wells$educ4)

fit_8 <- stan_glm(switch ~ c_dist100 + c_arsenic + c_educ4 +
                      c_dist100:c_educ4 + c_arsenic:c_educ4,
                  family = binomial(link="logit"), data = wells)

print(fit_8, digits=2)

stan_glm
 family:       binomial [logit]
 formula:      switch ~ 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.10 
c_arsenic          0.49   0.04 
c_educ4            0.19   0.04 
c_dist100:c_educ4  0.33   0.10 
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

LOO log score

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

Computed from 4000 by 3020 log-likelihood matrix

         Estimate   SE
elpd_loo  -1952.9 16.5
p_loo         6.6  0.3
looic      3905.8 33.0
------
Monte Carlo SE of elpd_loo is 0.0.

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.6       6.3  

loo_compare(loo7, loo8)

      elpd_diff se_diff
fit_8  0.0       0.0   
fit_7 -6.6       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$switch==1]-pred3[wells$switch==1],pred3[wells$switch==0]-pred8[wells$switch==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(switch ~ dist100 + log_arsenic, family = binomial(link = "logit"),
                   data = wells)

print(fit_3a, digits=2)

stan_glm
 family:       binomial [logit]
 formula:      switch ~ dist100 + log_arsenic
 observations: 3020
 predictors:   3
------
            Median MAD_SD
(Intercept)  0.52   0.06 
dist100     -0.98   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

LOO log score

(loo3a <- loo(fit_3a))

Computed from 4000 by 3020 log-likelihood matrix

         Estimate   SE
elpd_loo  -1952.2 16.3
p_loo         3.0  0.1
looic      3904.4 32.6
------
Monte Carlo SE of elpd_loo is 0.0.

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.3       4.4  

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

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

print(fit_4a, digits=2)

stan_glm
 family:       binomial [logit]
 formula:      switch ~ dist100 + log_arsenic + dist100:log_arsenic
 observations: 3020
 predictors:   4
------
                    Median MAD_SD
(Intercept)          0.49   0.07 
dist100             -0.88   0.14 
log_arsenic          0.99   0.11 
dist100:log_arsenic -0.23   0.18 

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

LOO log score

(loo4a <- loo(fit_4a))

Computed from 4000 by 3020 log-likelihood matrix

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

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.3   

Add interactions with education

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

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

print(fit_8a, digits=2)

stan_glm
 family:       binomial [logit]
 formula:      switch ~ 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.00   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.06   0.07 

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

LOO log score

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

Computed from 4000 by 3020 log-likelihood matrix

         Estimate   SE
elpd_loo  -1938.0 17.1
p_loo         6.1  0.2
looic      3875.9 34.3
------
Monte Carlo SE of elpd_loo is 0.0.

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.9       4.3