Fitting the same regression to many datasets. See Chapter 10 in Regression and Other Stories.


Load packages

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

Load data

data <- read.table(root("NES/data","nes.txt"))
head(data)
    year resid weight1 weight2 weight3 age gender race educ1 urban region
536 1952     1       1       1       1  25      2    1     2     2      1
537 1952     2       1       1       1  33      2    1     1     2      1
538 1952     3       1       1       1  26      2    1     2     2      1
539 1952     4       1       1       1  63      1    1     2     2      1
540 1952     5       1       1       1  66      2    1     2     2      2
541 1952     6       1       1       1  48      2    1     2     2      2
    income occup1 union religion educ2 educ3 martial_status occup2 icpsr_cty
536      4      2     1        1     3     3              1      2        NA
537      4      6     1        1     1     1              1      6        NA
538      3      6     2        2     3     3              1      6        NA
539      3      3     1        1     2     2              1      3        NA
540      1      6     2        1     4     4              1      6        NA
541      4      6     1        1     2     2              1      6        NA
    fips_cty partyid7 partyid3 partyid3_b str_partyid father_party mother_party
536       NA        6        3          3           3            3            3
537       NA        5        3          3           2            2            2
538       NA        4        2          2           1            1            1
539       NA        7        3          3           4            1           NA
540       NA        7        3          3           4            1            1
541       NA        3        1          1           2            1            1
    dlikes rlikes dem_therm rep_therm regis vote regisvote presvote
536      0      1        NA        NA     2    2         3        2
537     -1      3        NA        NA     2    2         3        1
538      0      5        NA        NA     2    2         3        2
539     -1      3        NA        NA     1    2         3        2
540     -2      0        NA        NA     2    2         3        2
541      0      4        NA        NA     2    2         3        2
    presvote_2party presvote_intent ideo_feel ideo7 ideo cd state inter_pre
536               2               2        NA    NA   NA NA    13        50
537               1               2        NA    NA   NA NA    13        50
538               2               2        NA    NA   NA NA    13        50
539               2               2        NA    NA   NA NA    13        50
540               2               2        NA    NA   NA NA    24        49
541               2               2        NA    NA   NA NA    24        49
    inter_post black female age_sq rep_presvote rep_pres_intent south real_ideo
536         NA     0      1    625            1               1     0        NA
537         NA     0      1   1089            0               1     0        NA
538         NA     0      1    676            1               1     0        NA
539         NA     0      0   3969            1               1     0        NA
540         NA     0      1   4356            1               1     0        NA
541         NA     0      1   2304            1               1     0        NA
    presapprov perfin1 perfin2 perfin presadm age_10 age_sq_10 newfathe newmoth
536         NA      NA      NA     NA      -1    2.5  6.250000        1       1
537         NA      NA      NA     NA      -1    3.3 10.889999        0       0
538         NA      NA      NA     NA      -1    2.6  6.759999       -1      -1
539         NA      NA      NA     NA      -1    6.3 39.690002       -1      NA
540         NA      NA      NA     NA      -1    6.6 43.559998       -1      -1
541         NA      NA      NA     NA      -1    4.8 23.040001       -1      -1
    parent_party white year_new income_new   age_new vote.1 age_discrete
536            2     1        1          1 -2.052455      1            1
537            0     1        1          1 -1.252455      1            2
538           -2     1        1          0 -1.952455      1            1
539           NA     1        1          0  1.747545      1            3
540           -2     1        1         -2  2.047545      1            4
541           -2     1        1          1  0.247545      1            3
    race_adj dvote rvote
536        1     0     1
537        1     1     0
538        1     0     1
539        1     0     1
540        1     0     1
541        1     0     1

Partyid model to illustrate repeated model use (secret weapon)

regress_year <- function (yr) {
  this_year <- data[data$year==yr,]
  fit <- stan_glm(partyid7 ~ real_ideo + race_adj + factor(age_discrete) +
                      educ1 + female + income,
                  data=this_year, warmup = 500, iter = 1500, refresh = 0,
                  save_warmup = FALSE, cores = 1, open_progress = FALSE)
  coefs <- cbind(coef(fit),se(fit))
}
summary <- array (NA, c(9,2,8))
for (yr in seq(1972,2000,4)){
  i <- (yr-1968)/4
  summary[,,i] <- regress_year(yr)
}
yrs <- seq(1972,2000,4)

Plot

coef_names <- c("Intercept", "Ideology", "Black", "Age_30_44", "Age_45_64", "Age_65_up", "Education", "Female", "Income")
par(mfrow=c(2,5), mar=c(2,3,2,2), tck=-.02, mgp=c(2,.7,0))
for (k in 1:9){
  plot(range(yrs), range(0,summary[k,1,]+.67*summary[k,2,],summary[k,1,]-.67*summary[k,2,]), type="n", xlab="", ylab="Coefficient", main=coef_names[k], mgp=c(1.2,.2,0), cex.main=1, cex.axis=1, cex.lab=1, tcl=-.1, bty="l", xaxt="n")
  axis(1, c(1972,1986,2000), mgp=c(.5,.3,0))
  abline(0,0, lty=2)
  points(yrs, summary[k,1,], pch=20)
  segments(yrs, summary[k,1,]-.67*summary[k,2,], yrs, summary[k,1,]+.67*summary[k,2,])
}

LS0tCnRpdGxlOiAiUmVncmVzc2lvbiBhbmQgT3RoZXIgU3RvcmllczogTmF0aW9uYWwgZWxlY3Rpb24gc3R1ZHkiCmF1dGhvcjogIkFuZHJldyBHZWxtYW4sIEplbm5pZmVyIEhpbGwsIEFraSBWZWh0YXJpIgpkYXRlOiAiYHIgZm9ybWF0KFN5cy5EYXRlKCkpYCIKb3V0cHV0OgogIGh0bWxfZG9jdW1lbnQ6CiAgICB0aGVtZTogcmVhZGFibGUKICAgIHRvYzogdHJ1ZQogICAgdG9jX2RlcHRoOiAyCiAgICB0b2NfZmxvYXQ6IHRydWUKICAgIGNvZGVfZG93bmxvYWQ6IHRydWUKLS0tCkZpdHRpbmcgdGhlIHNhbWUgcmVncmVzc2lvbiB0byBtYW55IGRhdGFzZXRzLiBTZWUgQ2hhcHRlciAxMCBpbgpSZWdyZXNzaW9uIGFuZCBPdGhlciBTdG9yaWVzLgoKLS0tLS0tLS0tLS0tLQoKCmBgYHtyIHNldHVwLCBpbmNsdWRlPUZBTFNFfQprbml0cjo6b3B0c19jaHVuayRzZXQobWVzc2FnZT1GQUxTRSwgZXJyb3I9RkFMU0UsIHdhcm5pbmc9RkFMU0UsIGNvbW1lbnQ9TkEpCiMgc3dpdGNoIHRoaXMgdG8gVFJVRSB0byBzYXZlIGZpZ3VyZXMgaW4gc2VwYXJhdGUgZmlsZXMKc2F2ZWZpZ3MgPC0gRkFMU0UKYGBgCgojIyMjIExvYWQgcGFja2FnZXMKCmBgYHtyIH0KbGlicmFyeSgicnByb2pyb290IikKcm9vdDwtaGFzX2ZpbGUoIi5ST1MtRXhhbXBsZXMtcm9vdCIpJG1ha2VfZml4X2ZpbGUoKQpsaWJyYXJ5KCJyc3RhbmFybSIpCmxpYnJhcnkoImZvcmVpZ24iKQpgYGAKCiMjIyMgTG9hZCBkYXRhCgpgYGB7ciB9CmRhdGEgPC0gcmVhZC50YWJsZShyb290KCJORVMvZGF0YSIsIm5lcy50eHQiKSkKaGVhZChkYXRhKQpgYGAKCiMjIyMgUGFydHlpZCBtb2RlbCB0byBpbGx1c3RyYXRlIHJlcGVhdGVkIG1vZGVsIHVzZSAoc2VjcmV0IHdlYXBvbikKCmBgYHtyIH0KcmVncmVzc195ZWFyIDwtIGZ1bmN0aW9uICh5cikgewogIHRoaXNfeWVhciA8LSBkYXRhW2RhdGEkeWVhcj09eXIsXQogIGZpdCA8LSBzdGFuX2dsbShwYXJ0eWlkNyB+IHJlYWxfaWRlbyArIHJhY2VfYWRqICsgZmFjdG9yKGFnZV9kaXNjcmV0ZSkgKwogICAgICAgICAgICAgICAgICAgICAgZWR1YzEgKyBmZW1hbGUgKyBpbmNvbWUsCiAgICAgICAgICAgICAgICAgIGRhdGE9dGhpc195ZWFyLCB3YXJtdXAgPSA1MDAsIGl0ZXIgPSAxNTAwLCByZWZyZXNoID0gMCwKICAgICAgICAgICAgICAgICAgc2F2ZV93YXJtdXAgPSBGQUxTRSwgY29yZXMgPSAxLCBvcGVuX3Byb2dyZXNzID0gRkFMU0UpCiAgY29lZnMgPC0gY2JpbmQoY29lZihmaXQpLHNlKGZpdCkpCn0Kc3VtbWFyeSA8LSBhcnJheSAoTkEsIGMoOSwyLDgpKQpmb3IgKHlyIGluIHNlcSgxOTcyLDIwMDAsNCkpewogIGkgPC0gKHlyLTE5NjgpLzQKICBzdW1tYXJ5WywsaV0gPC0gcmVncmVzc195ZWFyKHlyKQp9CnlycyA8LSBzZXEoMTk3MiwyMDAwLDQpCmBgYAoKIyMjIyBQbG90CgpgYGB7ciBldmFsPUZBTFNFLCBpbmNsdWRlPUZBTFNFfQppZiAoc2F2ZWZpZ3MpIHBkZihyb290KCJORVMvZmlncyIsInBhcnR5aWRfMS5wZGYiKSwgaGVpZ2h0PTIuNSwgd2lkdGg9Ny41KQpgYGAKYGBge3IgfQpjb2VmX25hbWVzIDwtIGMoIkludGVyY2VwdCIsICJJZGVvbG9neSIsICJCbGFjayIsICJBZ2VfMzBfNDQiLCAiQWdlXzQ1XzY0IiwgIkFnZV82NV91cCIsICJFZHVjYXRpb24iLCAiRmVtYWxlIiwgIkluY29tZSIpCnBhcihtZnJvdz1jKDIsNSksIG1hcj1jKDIsMywyLDIpLCB0Y2s9LS4wMiwgbWdwPWMoMiwuNywwKSkKZm9yIChrIGluIDE6OSl7CiAgcGxvdChyYW5nZSh5cnMpLCByYW5nZSgwLHN1bW1hcnlbaywxLF0rLjY3KnN1bW1hcnlbaywyLF0sc3VtbWFyeVtrLDEsXS0uNjcqc3VtbWFyeVtrLDIsXSksIHR5cGU9Im4iLCB4bGFiPSIiLCB5bGFiPSJDb2VmZmljaWVudCIsIG1haW49Y29lZl9uYW1lc1trXSwgbWdwPWMoMS4yLC4yLDApLCBjZXgubWFpbj0xLCBjZXguYXhpcz0xLCBjZXgubGFiPTEsIHRjbD0tLjEsIGJ0eT0ibCIsIHhheHQ9Im4iKQogIGF4aXMoMSwgYygxOTcyLDE5ODYsMjAwMCksIG1ncD1jKC41LC4zLDApKQogIGFibGluZSgwLDAsIGx0eT0yKQogIHBvaW50cyh5cnMsIHN1bW1hcnlbaywxLF0sIHBjaD0yMCkKICBzZWdtZW50cyh5cnMsIHN1bW1hcnlbaywxLF0tLjY3KnN1bW1hcnlbaywyLF0sIHlycywgc3VtbWFyeVtrLDEsXSsuNjcqc3VtbWFyeVtrLDIsXSkKfQpgYGAKYGBge3IgZXZhbD1GQUxTRSwgaW5jbHVkZT1GQUxTRX0KaWYgKHNhdmVmaWdzKSBkZXYub2ZmKCkKYGBgCgo=