Replication.Rmd

---
title: "Replication of Demographic Imperative Project"
author: "Grace Tien and Jessamin Birdsall"
date: "April 26, 2015"
output: html_document
---

The paper we have chosen to replicate is: "The Demographic Imperative in Religious Change in the United States," written by Hout, Greeley, and Wilde (2001). The question that the authors are trying to answer is: Over the course of the last century, why has affiliation with "mainline" Protestantism declined, while affiliation with "conservative" Protestantism has increased in the United States?  (Broadly speaking, "conservative" refers to more theologically conservative denominations and "mainline" refers to more liberal denominations.) To date, almost all sociological scholarship related to this question has assumed that mainline decline is the result of people converting/switching from mainline to conservative Protestantism. This paper argues that denominational switching actually accounts for *none* of the mainline decline. They argue that the most significant variable explaning mainline decline is *higher fertility* among conservative Protestants. In this paper, they use GSS data to test 4 hypotheses about the source of mainline decline: 1) higher fertility of conservative Protestants, 2) higher rate of mainline to conservative switching, 3) higher rate of apostasy among mainliners, 4) higher inflow of outsiders to conservative Protestant denominations. Looking at all cohorts born between 1900 and 1973, the authors develop demographic simulation models based on the variables of fertility, religious origins, and current religion. They then make counterfactual predictions in order to isolate the relative contributions of each of these variables to the observed decline in mainline affiliation. 

```{r}
# Install packages and load libraries. Because we are working with the GSS Stata file (downloaded from the NORC website), we have to use the foreign package in order to read it in R. To read the excel files, we need the xlsx package. 

library(dplyr)
library(ggplot2)
library(broom)
library(grid)
library(stargazer)
library(foreign)
library(xlsx)
library(gdata)

# Open the file. Within the code, we insert convert.underscore = TRUE because R Stata uses underscores in variable names, but R does not. 

GSS.df <- read.dta("GSS7214_R1.DTA", convert.factors = FALSE, convert.underscore=TRUE)

# Load the other documents shared by the authors of the paper. The authors shared a few excel files with us that show some of their calculations. We load them here for reference, although we are not able to interpret all of the calculations within them. 

maintrend.df <- read.dta("Hout/maintrend.dta", convert.factors = FALSE, convert.underscore=TRUE)

maincohort.df <- read.dta("Hout/maincohort.dta", convert.factors = F, convert.underscore=TRUE)

moblow.df <- read.dta("Hout/moblow.dta", convert.factors = F, convert.underscore=TRUE)

hout.df <- read.xlsx("Hout/Scenarios.xls", sheetName = "switching.within.prot")

# Select relevant variables. We are interested primarily in the variables of adult religious affiliation (relig), adult Protestant affiliation (fund), childhood Protestant affiliation (fund16), and total number of children born (childs) to the women in the sample. We filter out data collected after the year 1998, because the authors the used GSS data until 1998. Like the authors, we consider only women aged 25-74 who were born between 1903 and 1973.

GSS1.df <- GSS.df %>%
  select(cohort, age, sex, relig, relig16, fund, fund16, childs, year) %>%
  arrange(desc(cohort)) %>%
  filter(year <= 1998, cohort >= 1903, cohort <= 1973, age >= 25, age <= 74, !is.na(relig), !is.na(relig16), !is.na(fund), !is.na(fund16), !is.na(childs))
```

*Figure 1: Fertility by Denominational Type* 

```{r}
# Replicate Figure 1. The first figure in the paper is a graph of fertility by denominational type (e.g. mainline and conservative). The x axis is year of birth (ranging from 1903 to 1950) and the y axis is births per woman.

# Recode denominations. Before we can do the calculations that go into this figure, we have to first recode the denominations from the GSS variables into values we can use. A person is considered mainline if s/he answered "Protestant" to the relig question and either "liberal" or "moderate" to the fund question. This is consistent with the authors' categorization, as best we can understand it. 

GSS1.df$adult_mainline <- ifelse((GSS1.df$fund == 2 | GSS1.df$fund == 3) & GSS1.df$relig == 1, 1, 0)
adult_mainline <- GSS1.df$adult_mainline

GSS1.df$adult_cons <- ifelse((GSS1.df$fund == 1) & GSS1.df$relig == 1, 1, 0)
adult_cons <- GSS1.df$adult_cons

GSS1.df$adult_prot <- ifelse((GSS1.df$fund <=3 ) & GSS1.df$relig == 1, 1, 0)
adult_prot <- GSS1.df$adult_prot

# Subset the mainline fertility data.

mainline_childs <- GSS1.df  %>% 
  select(cohort, childs)  %>% 
  filter(adult_mainline == 1) %>%
  group_by(cohort) %>%
  summarise(mean(childs)/n_distinct(cohort))

# Subset the conservative fertility data.

cons_childs <- GSS1.df  %>% 
  select(cohort, childs)  %>% 
  filter(adult_cons == 1) %>%
  group_by(cohort) %>%
  summarise(mean(childs)/n_distinct(cohort))

# Combined fertility rates for Protestants.

all_childs <- GSS1.df  %>% 
  select(cohort, childs)  %>% 
  filter(adult_prot == 1) %>%
  group_by(cohort) %>%
  summarise(mean(childs)/n_distinct(cohort))

# Convert to dataframes that ggplot can recognize.

# Children born to mainline women. "CEB" stands for "children ever born."

m_childs <- as.data.frame(mainline_childs)
names(m_childs) = c("m_cohort","m_CEB")

# Children born to conservative women:

c_childs <- as.data.frame(cons_childs)
names(c_childs) = c("c_cohort", "c_CEB")

# Children born to all Protestant women:

all_childs <- as.data.frame(all_childs)
names(all_childs) = c("a_cohort", "a_CEB")

# Plot CEB (children ever born) for all 3 categories of women (mainline, conservative, and all together) by year of birth. The authors use "loess" rather than linear regression to fit the line more closely to the data. Shape = 1 makes the dots empty rather than filled in. 

p1 <- ggplot() +
  geom_point(data =  m_childs, aes(m_cohort, m_CEB), color = "blue") + 
  geom_smooth(data = m_childs, aes(m_cohort, m_CEB), size = 1.5, method = "loess", span = .2, se = FALSE, color = "blue") +
  ggtitle("Fertility by Denomination Type") + 
  labs(x = "Year of Birth", y = "Births per Woman") +
  ylim(c(1, 4)) + 
  xlim(c(1903, 1950)) +
  geom_point(data = c_childs, aes(c_cohort, c_CEB), shape = 1, color = "red") + 
  geom_smooth(data = c_childs, aes(c_cohort, c_CEB), method = "loess", span = .2, se = FALSE, color = "red") +
  geom_smooth(data = all_childs, aes(a_cohort, a_CEB), method = "loess", span = .2, se = FALSE, color = "black")
# annotation_custom(grobTree(textGrob("FIG 1 - Fertility (children ever born) by denomination type: US Protestant women, 45-69 years old. CEB data are smoothed using loess regression.", x = 0.1, y = 0.1, hjust = 0, gp = gpar(col = "grey60", fontsize = 10, fontface = "italic")))) 
# For some reason this annotation is not showing up properly. 
# The red line corresponds to mainliners, blue to conservatives, and black to the combined group. 
p1

```

*Figure 2: Proportion of Protestants professing a mainline denomination by year*

```{r}
# Replicate Figure 2. This figure shows the proportion of U.S. Protestants who identify as mainline by year in which the survey was taken. The x axis is year of survey (ranging from 1973 to 1998), and the y axis is proportion of Protestants professing a mainline affiliation. 

# Calculate the proportion of Mainline Protestants by year.

m_prop <- GSS1.df  %>% 
  group_by(year) %>%
  summarise(total.prot = sum(adult_prot), total.mainline= sum(adult_mainline)) %>%
  mutate(prop.m = total.mainline/total.prot)

# Calculate confidence intervals

m_prop$se <- sqrt((m_prop$prop.m * (1-m_prop$prop.m))/m_prop$total.prot)
m_prop$lower <- m_prop$prop.m - 1.96 * m_prop$se
m_prop$upper <- m_prop$prop.m + 1.96 * m_prop$se

# Plot.

p2 <- ggplot(data = m_prop, aes(x = year, y = prop.m)) + geom_point(color = "firebrick") + ggtitle("Proportion of Mainline Protestants") +
  labs(x = "Year", y = "Percentage Mainline") +
  ylim(c(0.3, 0.6)) +
  xlim(c(1970, 2000)) +
  geom_smooth(method = "loess", span = .8, se = F) +
  geom_errorbar(aes(x = year, ymin = lower, ymax = upper), colour = "black")
p2
  
```

*Figure 3: Proportion of Protestants Professing a Mainline Denomination of Year of Birth*

```{r}
# This figure shows the proportion of U.S. Protestants who identify as mainline by cohort. The x axis is year of birth (ranging from 1903 to 1973), and the y axis is proportion of Protestants professing a mainline affiliation. 

# Calculate the proportion of mainline Protestants out of the total Protestants. Group by cohort. 

m_prop_cohort <- GSS1.df  %>% 
  group_by(cohort) %>%
  summarise(total.prot = sum(adult_prot), total.mainline= sum(adult_mainline)) %>%
  mutate(prop.m = total.mainline/total.prot)

# Calculate confidence intervals

m_prop_cohort$se <- sqrt((m_prop_cohort$prop.m * (1-m_prop_cohort$prop.m))/m_prop_cohort$total.prot)
m_prop_cohort$lower <- m_prop_cohort$prop.m - 1.96 * m_prop_cohort$se
m_prop_cohort$upper <- m_prop_cohort$prop.m + 1.96 * m_prop_cohort$se

# Plot

p3 <- ggplot(data = m_prop_cohort, aes(x = cohort, y = prop.m)) + geom_point(color = "firebrick") + ggtitle("Proportion of Mainline Protestants by Year of Birth") +
  labs(x = "Year of Birth", y = "Percentage Mainline") +
  ylim(c(0.2, 0.8)) +
  xlim(c(1900, 1975)) +
  geom_smooth(method = "lm", se = F) +
  geom_smooth(method = "loess", span = .4, se = F) +
  geom_errorbar(aes(x = cohort, ymin = lower, ymax = upper), colour = "black")
p3

```

*Figure 5: Observed Percentage of 25-74-year-old U.S. Protestants professing a mainline denomination and that predicted by switching between Protestant denominations by year of birth*

```{r}
# The purpose of this figure is to show the observed decline in mainline affiliation, and the predicted decline if we look only at switching. We draw upon the data stored in the hout.df because it is otherwise unclear how the authors generated switching from the raw data. 

# Model 1: this will help us to predict the proportion of mainline affiliation by proportion of people who have switched from mainline to conservative. 

model1 <- glm(prop.main ~ prop.switch.mtoc, data = hout.df)
stargazer(model1, type = "text")

# Model 2: this will help us to predict the proportion of mainline affiliation by proportion of people who have switched from conservative to mainline.

model2 <- glm(prop.main ~ prop.switch.ctom, data = hout.df)
stargazer(model2, type = "text")

# Calculate fitted values for Model 1

pred_prop1 <- predict(model1, type = "response", newdata = hout.df)
hout.df$pred_prop1 <- pred_prop1

# Calculate fitted values for Model 2

pred_prop2 <- predict(model2, type = "response", newdata = hout.df)
hout.df$pred_prop2 <- pred_prop2

# Plot observed and fitted values. Smooth using bandwidth of .4

ggplot(hout.df, aes(x = birth.year, y = prop.main)) +
  geom_point(color = "black") +
  geom_smooth(aes(y = prop.main), method = "loess", se = F, span = .4, color = "black") +
  geom_point(aes(y = pred_prop1), color = "red") +
  geom_smooth(aes(y = pred_prop1), method = "loess", se = F, span = .4, color = "red") +
  geom_point(aes(y = pred_prop2), color = "blue") +
  geom_smooth(aes(y = pred_prop2), method = "loess", se = F, span = .4, color = "blue") +
  ylim(.3, .7) +
  ggtitle("Proportion Mainline, Observed and Predicted by Switching") +
  labs(x = "Year of Birth", y = "Percentage Mainline")

# The black line corresponds to the observed decline, the red line to the predictions based on Model 1 (mainline to conservative switching), and the blue line to the predictions based on Model 2 (conservative to mainline switching). Model 1 is producing predictions in line with the authors' observations, but Model 2 is not. We structured the model in the same way for both, so we are not sure why Model 2 is not producing the same predictions. 

```

*Table 1: Denominational Switching by Type of Origin Denomination*

```{r}
# As mentioned above, many social scientists and religious leaders have assumed that the phenonemon of switching (converting from one denomination to the other) is the primary explanatory variable behind the decline in mainline Protestant affiliation. The purpose of this table is to capture to what extent switching has happened across time, both within Protestantism and outside (to other religions or to no religion).

# Create new data frame with mutated variables that capture switching dynamics. 

GSS2.df <- GSS1.df %>% 
  select(cohort, relig, relig16, fund, fund16) %>%
  mutate(adult_mainline = ifelse(fund == 2 | fund == 3 & relig == 1, 1, 0), 
kid_mainline = ifelse(fund16 == 2 | fund16 == 3 & relig16 == 1, 1, 0), 
adult_other = ifelse(fund >2 & relig == 2 | relig == 3 | relig == 5 | relig == 6 | relig == 7 | relig == 8 | relig == 9 | relig == 10 | relig == 11 | relig == 12 | relig == 13, 1, 0),
adult_none = ifelse(fund >1 & relig == 4, 1, 0),
pm_same = sum(ifelse(adult_mainline==1 & kid_mainline == 1, 1, 0))/sum(ifelse(kid_mainline==1,1,0)),
         pm_otherprot = sum(ifelse(adult_mainline==0 & kid_mainline == 1,1,0))/sum(ifelse(kid_mainline==1,1,0)),
         pm_otherrel = sum(ifelse(adult_other==1 & kid_mainline == 1,1,0))/sum(ifelse(kid_mainline==1,1,0)),
         pm_none = sum(ifelse(adult_none==1 & kid_mainline == 1,1,0))/sum(ifelse(kid_mainline==1,1,0)),
         pc_same = sum(ifelse(adult_mainline==0 & kid_mainline == 0,1,0))/sum(ifelse(kid_mainline==0,1,0)),
         pc_otherprot = sum(ifelse(adult_mainline==1 & kid_mainline == 0,1,0))/sum(ifelse(kid_mainline==0,1,0)),
         pc_otherrel = sum(ifelse(adult_other==1 & kid_mainline == 0,1,0))/sum(ifelse(kid_mainline==0,1,0)),
         pc_none = sum(ifelse(adult_none==1 & kid_mainline == 0,1,0))/sum(ifelse(kid_mainline==0,1,0)))

# Put in table format that aggregates into decades.

tb.1 <- GSS2.df %>%
  mutate(tb_period = ifelse(cohort >= 1900 & cohort <= 1909, "1900-9", "GSS2.df"),
       tb_period = ifelse(cohort >= 1910 & cohort <= 1919, "1910-19", tb_period),
       tb_period = ifelse(cohort >= 1920 & cohort <= 1929, "1920-29", tb_period),
       tb_period = ifelse(cohort >= 1930 & cohort <= 1939, "1930-39", tb_period),
       tb_period = ifelse(cohort >= 1940 & cohort <= 1940, "1940-49", tb_period),
       tb_period = ifelse(cohort >= 1950 & cohort <= 1959, "1950-59", tb_period),
       tb_period = ifelse(cohort >= 1960 & cohort <= 1973, "1960-73", tb_period)) %>%
  arrange(cohort) %>% 
  group_by(tb_period) %>%
  summarise(mean_pm_same = mean(pm_same),
            mean_pm_otherprot = mean(pm_otherprot),
            mean_pm_otherrel = mean(pm_otherrel),
            mean_pm_none = mean(pm_none),
            mean_pc_same = mean(pc_same),
            mean_pc_otherprot = mean(pc_otherprot),
            mean_pc_otherrel = mean(pc_otherrel),
            mean_pc_none = mean(pc_none)) 

tb.1

# After many, many hours, this code is still not doing exactly what we want it to do. In an earlier version, we were able to get the code to calculate the proportions by decade, but for some reason the code is now just returning one proporation across all decades. 
```

*Table 2: Denominational Switching by Type of Current Denomination*

Table 2 is similar to Table 1, but in reverse. It shows, for people who currently identify as mainline, the proportions who identified as other Protestant, other religion, or no religion as a child. It then shows the same for people who currently identify as conservative.

```{r, eval = FALSE}
# We have asked R not to evaluate this chunk of code, because, after repeated attempts to knit, we kept getting the following error message: Error in parse(text = x, srcfile = src) : <text>:22:1:unexpected symbol 21: 22: tb.2 ^ Calls: <Anonymous>... evaluate -> parse_all -> parse_all.character -> parse

# Create new data frame with mutated variables that capture switching dynamics
library(dplyr)
GSS3.df <- GSS1.df %>% 
  select(cohort, relig, relig16, fund, fund16) %>%
  mutate(adult_mainline = ifelse(fund == 2 | fund == 3 & relig == 1, 1, 0), 
  adult_cons = ifelse(fund == 1 & relig == 1, 1, 0),
  kid_cons = ifelse(fund16 == 1 & relig16 == 1, 1, 0),
kid_mainline = ifelse(fund16 == 2 | fund16 == 3 & relig16 == 1, 1, 0), 
kid_other = ifelse(fund16 > 2 & relig16 == 2 | relig16 == 3 | relig16 == 5 | relig16 == 6 | relig16 == 7 | relig16 == 8 | relig16 == 9 | relig16 == 10 | relig16 == 11 | relig16 == 12 | relig16 == 13, 1, 0),
kid_none = ifelse(fund16 >1 & relig16 == 4, 1, 0),
pm_same = sum(ifelse(kid_mainline == 1 & adult_mainline==1, 1, 0))/sum(ifelse(adult_mainline == 1,1,0)),
         pm_otherprot = sum(ifelse(kid_cons == 1 & adult_mainline==1,1,0))/sum(ifelse(adult_mainline == 1,1,0)),
         pm_otherrel = sum(ifelse(kid_other == 1 & adult_mainline == 1,1,0))/sum(ifelse(adult_mainline == 1,1,0)),
         pm_none = sum(ifelse(kid_none ==1 & adult_mainline == 1,1,0))/sum(ifelse(adult_mainline == 1,1,0)),
         pc_same = sum(ifelse(kid_cons == 1 & adult_cons == 1,1,0))/sum(ifelse(adult_cons == 1,1,0)),
         pc_otherprot = sum(ifelse(kid_cons == 1& adult_cons == 1,1,0))/sum(ifelse(adult_cons == 1,1,0)),
         pc_otherrel = sum(ifelse(kid_other == 1 & adult_cons == 1,1,0)/sum(ifelse(adult_cons == 1,1,0)),
         pc_none = sum(ifelse(kid_none == 1 & adult_cons == 1,1,0)/sum(ifelse(adult_cons == 1,1,0)))

# Put in table format that aggregates into decades

tb.2 <- GSS3.df %>%
 arrange(cohort)  %>% 
  mutate(tb_period = ifelse(cohort >= 1900 & cohort <= 1909, "1900-9", "GSS3.df"),
       tb_period = ifelse(cohort >= 1910 & cohort <= 1919, "1910-19", tb_period),
       tb_period = ifelse(cohort >= 1920 & cohort <= 1929, "1920-29", tb_period),
       tb_period = ifelse(cohort >= 1930 & cohort <= 1939, "1930-39", tb_period),
       tb_period = ifelse(cohort >= 1940 & cohort <= 1940, "1940-49", tb_period),
       tb_period = ifelse(cohort >= 1950 & cohort <= 1959, "1950-59", tb_period),
       tb_period = ifelse(cohort >= 1960 & cohort <= 1973, "1960-73", tb_period)) %>%
  group_by(tb_period) %>%
  summarise(mean_pm_same = mean(pm_same),
            mean_pm_otherprot = mean(pm_otherprot),
            mean_pm_otherrel = mean(pm_otherrel),
            mean_pm_none = mean(pm_none),
            mean_pc_same = mean(pc_same),
            mean_pc_otherprot = mean(pc_otherprot),
            mean_pc_otherrel = mean(pc_otherrel),
            mean_pc_none = mean(pc_none)) 
#As above, this code is still not doing exactly what we want it to do. In an earlier version, we were #able to get the code to calculate the proportions by decade, but for some reason the code is now just #returning one proporation across all decades. 
#
tb.2

```

*Extension of Figure 1: Fertility by Denominational Type*
This paper was published in 2001, and so the authors used GSS data up through 1998 to analyze trends in denominational affiliation and fertility. They looked at cohorts born between 1903 and 1973. We have the the GSS cumulative data through 2006, so we are curious to see how, if at all, trends in affiliation or fertility may have changed since 1998. 

```{r}
# Extend the dataframe to include all data until 2006, and cohorts through 1981. 
extended.df <- GSS.df %>%
  select(cohort, age, sex, relig, relig16, fund, fund16, childs, year) %>%
  arrange(desc(cohort)) %>%
  filter(cohort >= 1903, cohort <= 1981, age >= 25, age <= 74, !is.na(relig), !is.na(relig16), !is.na(fund), !is.na(fund16), !is.na(childs), !is.na(cohort))

# Recode denominations. 

extended.df$adult_mainline2 <- ifelse((extended.df$fund == 2 | extended.df$fund == 3) & extended.df$relig == 1, 1, 0)
adult_mainline2 <- extended.df$adult_mainline2

# Subset the mainline fertility data.

mainline_childs2 <- extended.df  %>% 
  select(cohort, childs)  %>% 
  filter(adult_mainline2 == 1) %>%
  group_by(cohort) %>%
  summarise(mean(childs)/n_distinct(cohort))

# Subset the conservative fertility data.

cons_childs2 <- extended.df  %>% 
  select(cohort, childs)  %>% 
  filter(adult_mainline2 == 0) %>%
  group_by(cohort) %>%
  summarise(mean(childs)/n_distinct(cohort))

# Convert to dataframes that ggplot can recognize.

# Children born to mainline women. "CEB" stands for "children ever born."

m_childs2 <- as.data.frame(mainline_childs2)
names(m_childs2) = c("m_cohort","m_CEB")

# Children born to conservative women:

c_childs2 <- as.data.frame(cons_childs2)
names(c_childs2) = c("c_cohort", "c_CEB")

# Plot 

p.ext <- ggplot() +
  geom_point(data =  m_childs2, aes(m_cohort, m_CEB), color = "blue") + 
  geom_smooth(data = m_childs2, aes(m_cohort, m_CEB), size = 1.5, method = "loess", span = .2, se = FALSE, color = "blue") +
  ggtitle("Fertility by Denomination Type") + 
  labs(x = "Year of Birth", y = "Births per Woman") +
  ylim(c(1, 4)) + 
  xlim(c(1903, 1983)) +
  geom_point(data = c_childs2, aes(c_cohort, c_CEB), shape = 1, color = "red") + 
  geom_smooth(data = c_childs2, aes(c_cohort, c_CEB), method = "loess", span = .2, se = FALSE, color = "red") 
# annotation_custom(grobTree(textGrob("Extension Figure - Fertility (children ever born) by denomination type: US Protestant women, 45-69 years old. CEB data are smoothed using loess regression.", x = 0.1, y = 0.1, hjust = 0, gp = gpar(col = "grey60", fontsize = 10, fontface = "italic")))) 
# This annotation is not showing up properly. 
# The red line corresponds to mainliners and the blue to conservatives.  
p.ext

```
When we add in 8 more years of data and 8 more cohorts, we see that fertility has continued to decline for both mainline and conservative Protestants. Between 1973 and 1980, the slope of decline for conservatives appears to be steeper than that of mainlines, which suggests that perhaps conservative fertility rates may intersect with mainline fertility rates within the next decade. This would be an interesting phenomenon to investigate further in the future, as we have several hypotheses about why fertility rates between the two denominations may be converging (e.g. rising levels of education among conservatives, later age of marriage for conservative women).

sum <- 0

for(i in  1:50){

  sum <- sum + i

}