Regression and Other Stories: Gay

Simple models (linear and discretized age) and attitudes as a function of age. See Chapter 12 in Regression and Other Stories.

---

Load packages

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

Load data

data <- read.csv(root("Gay/data","naes04.csv"))
age <- seq(18,91)
response <- as.character(data[,"gayFavorStateMarriage"])
n_age <- length(age)
yes <- rep(NA, n_age)
no <- rep(NA, n_age)
for (i in 1:n_age){
  if (i == n_age)
    ok <- data[,"age"] >= age[i]
  else
    ok <- data[,"age"] == age[i]
  yes[i] <- sum(ok & response=="Yes", na.rm=TRUE)
  no[i] <- sum(ok & response=="No", na.rm=TRUE)
}

support <- yes/(yes + no)
age_cutpoints <- c(0, seq(29, 79, 10), 100)
age_discrete <- cut(age, age_cutpoints)
gay_total <- data.frame(support, age, age_discrete)

Fit linear regression model

fit_linear <- stan_glm(support ~ age, data=gay_total, refresh=0)
print(fit_linear)

stan_glm
 family:       gaussian [identity]
 formula:      support ~ age
 observations: 74
 predictors:   2
------
            Median MAD_SD
(Intercept) 0.6    0.0   
age         0.0    0.0   

Auxiliary parameter(s):
      Median MAD_SD
sigma 0.0    0.0   

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

Fit discretized age model

fit_binned <- stan_glm(support ~ factor(age_discrete), data=gay_total, refresh=0)
print(fit_binned)

stan_glm
 family:       gaussian [identity]
 formula:      support ~ factor(age_discrete)
 observations: 74
 predictors:   7
------
                             Median MAD_SD
(Intercept)                   0.5    0.0  
factor(age_discrete)(29,39]  -0.1    0.0  
factor(age_discrete)(39,49]  -0.1    0.0  
factor(age_discrete)(49,59]  -0.1    0.0  
factor(age_discrete)(59,69]  -0.2    0.0  
factor(age_discrete)(69,79]  -0.3    0.0  
factor(age_discrete)(79,100] -0.3    0.0  

Auxiliary parameter(s):
      Median MAD_SD
sigma 0.0    0.0   

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

Plot linear regression model

par(mar=c(3,3,01,0), mgp=c(1.7, .5, 0), tck=-.02)
plot(age, support, xlim=c(min(age)-1, max(age)+1), ylim=c(0, 0.65), xaxs="i", yaxs="i", xlab="Age", ylab="Support for same-sex marriage", xaxt="n", yaxt="n", bty="l", main="Linear regression", cex.main=1, type="n")
points(age, support, pch=20, col="gray40")
axis(1, seq(20, 80, 20))
axis(2, seq(0, 0.60, 0.20), paste(seq(0, 60, 20), "%", sep=""))
abline(coef(fit_linear)[1], coef(fit_linear)[2], lwd=2)

Plot discretized age model

par(mar=c(3,3,1,0), mgp=c(1.7, .5, 0), tck=-.02)
plot(age, support, xlim=c(min(age)-1, max(age)+1), ylim=c(0, 0.65), xaxs="i", yaxs="i", xlab="Age", ylab="Support for same-sex marriage", xaxt="n", yaxt="n", bty="l", main="Discretized age predictors", cex.main=1, type="n")
points(age, support, pch=20, col="gray40")
axis(1, seq(20, 80, 20))
axis(2, seq(0, 0.60, 0.20), paste(seq(0, 60, 20), "%", sep=""))
for (k in 1:(length(age_cutpoints) - 1)){
  lines(age_cutpoints[k:(k+1)], rep(coef(fit_binned)[1] + ifelse (k==1, 0, coef(fit_binned)[k]), 2), lwd=2)
}