#v Dijk 24-4-2025

####Preparation ------
library("eha") 
library("lubridate")
library("dplyr")
library("coxme")
library("reshape")
library("Hmisc")
library("tidyverse")
library("ggplot2")
library("gridExtra")
library("ggpubr")

getwd()
#read basic datafiles
#cleaned civil registry for Zeeland 2017 and file with all occupations: https://datasets.iisg.amsterdam/dataset.xhtml?persistentId=hdl:10622/QUJNCD
rm(list = ls())
setwd("HERE") #set working directory
pedigree <-  read.csv ("Pedigree.txt", header=T, sep=",",stringsAsFactors=F)
hisclass <- read.csv("Occupations.txt", sep =",", header=T)
pedigree[, c(11,23,24,25,32,39,44,49,54,59)] <- lapply(pedigree[, c(11,23,24,25,32,39,44,49,54,59)], as.Date, origin="1900-1-1")
str(pedigree)
head(pedigree$LastEntryDate,50)
head(pedigree$LastEntryDate,50)
matmort <-  read.csv("Matmortzeeland2.csv", header=T, sep=";",stringsAsFactors=F)

setwd("HERE") #set working directory

#check the birth years
table(year(pedigree$B_date))
#check the marriage years
table(year(pedigree$M_date_1))



#basic data construction-------------

#create basic dataframe with id, mother, father, sex
persons <- pedigree[,c(1:3,9)] 
colnames(persons) <- c("personid", "mother", "father", "sex")
length(which(!is.na(persons$mother) & is.na(persons$father))) #27,129

#add birthyear and deathyear
str(pedigree)
bdyear <- pedigree[,c("B_date", "D_date", "Id_person")]
colnames(bdyear) <- c("bdate", "ddate", "personid")
persons <- merge(persons, bdyear, all.x=T, by="personid")
str(persons)

#age at death
agedsg1 <- persons[,c("personid", "ddate", "bdate")]
agedsg1$aged <- as.numeric(difftime(agedsg1$ddate,agedsg1$bdate,units="weeks")) / 52
agedsg1 <- agedsg1[,c("personid", "aged")]
persons <- merge(persons, agedsg1, by="personid", all.x=T)
#NA in age at death: make it 999
length(which(is.na(persons$aged)))
persons$aged[is.na(persons$aged)] <- 999

#birthyear and deathyear
persons$byr <- (as.numeric(format(as.Date(persons$bdate, origin='1900-1-1'), '%Y')))
persons$dyr <- (as.numeric(format(as.Date(persons$ddate, origin='1900-1-1'), '%Y')))
barplot(table(persons$byr))
barplot(table(persons$dyr))

#last date of observation
lastentry <- pedigree[,c("LastEntryDate", "Id_person")]
colnames(lastentry) <- c("lldate", "personid")
persons <- merge(persons, lastentry, all.x=T, by="personid")
colnames(persons)

#mother's and father's bdate and ddate
agedg1 <- persons[,c("personid", "bdate", "lldate", "ddate")]
colnames(agedg1)[c(2,3,4)] <- c("bdatem", "lldatem", "ddatem")
persons <- merge(persons, agedg1, by.x="mother", by.y="personid", all.x=T)
colnames(agedg1)[c(2,3,4)] <- c("bdatef", "lldatef", "ddatef")
persons <- merge(persons, agedg1, by.x="father", by.y="personid", all.x=T)
str(persons)
length(which(is.na(persons$bdatem))) 
length(which(is.na(persons$lldatem))) 
length(which(is.na(persons$bdatef)))
length(which(is.na(persons$lldatef))) 

#time between birth and mother's death
str(persons)
persons$deathgapmother <- difftime(persons$ddatem,persons$bdate,units="days")
persons$deathgapmother[is.na(persons$deathgapmother)] <- 999999
table(persons$deathgapmother)
#maternal mortality: 0/1 variable
persons$matmort <- 0
persons$matmort[persons$deathgapmother < 43] <- 1
table(persons$matmort)
#paternal mortality
persons$deathgapfather <- difftime(persons$ddatef,persons$bdate,units="days")
persons$deathgapfather[is.na(persons$deathgapfather)] <- 999999
table(persons$deathgapfather)
persons$fatmort <- 0
persons$fatmort[persons$deathgapfather < 43] <- 1
table(persons$fatmort)

gc()




####Socioeconomic variables####
#SES is uncategorized. lower hisclass = higher occupational status. find highest
hisclass$Id_person <- as.numeric(hisclass$Id_person)
hisclass$HISCLASS <- as.numeric(hisclass$HISCLASS)
hisclasssort <- hisclass[with(hisclass, order(hisclass$Id_person, hisclass$HISCLASS)),]
hisclasssort <- hisclasssort[-which(duplicated(hisclasssort$Id_person)),]
hisclasssort <- hisclasssort[c("Id_person", "HISCLASS")]

#recode HISCLASS to HISCLASS-6 plus missing
table(hisclass$HISCLASS)
hisclasssort$his5 <- "0. Missing"
hisclasssort$his5[hisclasssort$HISCLASS==1] <- "1. Elite"
hisclasssort$his5[hisclasssort$HISCLASS==2] <- "1. Elite"
hisclasssort$his5[hisclasssort$HISCLASS==3] <- "2. Lower middle class"
hisclasssort$his5[hisclasssort$HISCLASS==4] <- "2. Lower middle class"
hisclasssort$his5[hisclasssort$HISCLASS==5] <- "2. Lower middle class"
hisclasssort$his5[hisclasssort$HISCLASS==6] <- "2. Lower middle class"
hisclasssort$his5[hisclasssort$HISCLASS==7] <- "3. Skilled workers"
hisclasssort$his5[hisclasssort$HISCLASS==8] <- "4. Farmers and fishermen"
hisclasssort$his5[hisclasssort$HISCLASS==9] <- "3. Skilled workers"
hisclasssort$his5[hisclasssort$HISCLASS==10] <- "5. Unskilled workers"
hisclasssort$his5[hisclasssort$HISCLASS==11] <- "5. Unskilled workers"
hisclasssort$his5[hisclasssort$HISCLASS==12] <- "5. Unskilled workers"
hisclasssort$his5[hisclasssort$HISCLASS==13] <- "5. Unskilled workers"
table(hisclasssort$his5)
hisclasssort <- hisclasssort[c("Id_person", "his5")]

#Merge hisclass of the father to the births
str(persons)
persons <- merge(persons, hisclasssort, by.x="father", by.y="Id_person", all.x=T)
persons$his5[is.na(persons$his5)] <- "0. Missing"
table(persons$his5)
str(persons)


#select the research group
#kids with known parent
anl <- subset(persons, !is.na(persons$mother))
#known follow up of mother after the birth of the child
str(anl)
anl$followup <- difftime(anl$lldatem, anl$bdate,units="days")
length(which(anl$followup > 0)) #13078 women with followup after birth, only?
#remove women that are not observed after birth but without dropping the women that die on the day of birth of children
anl <- subset(anl, anl$followup > 0 | anl$ddatem==anl$bdate)


#kids born after 1812 and before 1913
anl <- subset(anl, anl$byr > 1811)
anl <- subset(anl, anl$byr < 1914)

#women aged 15 to 45 (at the time of birth)
anl$bagem <- difftime(anl$bdate, anl$bdatem,units="days") / 365.25
anl <- subset(anl, anl$bagem > 15)
anl <- subset(anl, anl$bagem < 46)

#by period: 
anl$bcohort <- 0
anl$bcohort[anl$byr > 1811 & anl$byr < 1840] <- 1
anl$bcohort[anl$byr > 1839 & anl$byr < 1870] <- 2
anl$bcohort[anl$byr > 1869 & anl$byr < 1890] <- 3
anl$bcohort[anl$byr > 1889 & anl$byr < 1914] <- 4
table(anl$bcohort)
table(anl$bcohort, anl$matmort)

#analyse
anl$start <- 0
anl$end <- 42
anl$end[anl$deathgapmother < 42] <- anl$deathgapmother[anl$deathgapmother < 42]
anl$event <- anl$matmort

#numbers
str(anl)
#mothers
length(unique(anl$mother))
#children
length(unique(anl$personid))
#events
table(anl$event)
table(anl$bcohort, anl$event)

#check for missings
length(which(is.na(anl$sex)))
length(which(is.na(anl$his5)))
length(which(is.na(anl$byr)))
length(which(is.na(anl$start)))
length(which(is.na(anl$end)))
length(which(is.na(anl$event)))
length(which(anl$start==anl$end))
table(anl$end)
anl$end[anl$start==anl$end] <- 0.25

#set reference category
anl$his5 <- factor(anl$his5, levels=c("5. Unskilled workers", "1. Elite", "2. Lower middle class", "3. Skilled workers", "4. Farmers and fishermen", "0. Missing"))

#make a separate dataset with observations on mothers, not children
vecmothers <- unique(anl$mother)
anl2 <- subset(persons, persons$personid %in% vecmothers)

#check for how many mothers there is a date of death or an observation after age 45
length(which(!is.na(anl2$ddate)))
length(which(!is.na(anl2$ddate))) / length(unique(anl2$personid))

#follow-up observation
anl2$llage <- difftime(anl2$lldate, anl2$bdate,units="days") / 365.25
length(which(!is.na(anl2$ddate) | anl2$llage > 45 ))
length(which(!is.na(anl2$ddate) | anl2$llage > 45 )) / length(unique(anl2$personid))

#analyse: add rows for each of the 4 periods 
anl2$bcohort <- 0

#extra data frames
extrarows1 <- anl2
extrarows2 <- anl2
extrarows3 <- anl2
extrarows4 <- anl2

#define birth cohort
extrarows1$bcohort <- 1
extrarows2$bcohort <- 2
extrarows3$bcohort <- 3
extrarows4$bcohort <- 4

#starting age in each of these periods
extrarows1$startage <- as.numeric(difftime(as.Date("1812-01-01"),extrarows1$bdate,units="days") / 365.25)
extrarows2$startage <- as.numeric(difftime(as.Date("1840-01-01"),extrarows2$bdate,units="days") / 365.25)
extrarows3$startage <- as.numeric(difftime(as.Date("1870-01-01"),extrarows3$bdate,units="days") / 365.25)
extrarows4$startage <- as.numeric(difftime(as.Date("1890-01-01"),extrarows4$bdate,units="days") / 365.25)

#if start age is below 20, censor at age 20
extrarows1$startage[extrarows1$startage < 20] <- 20
extrarows2$startage[extrarows2$startage < 20] <- 20
extrarows3$startage[extrarows3$startage < 20] <- 20
extrarows4$startage[extrarows4$startage < 20] <- 20

#ending age at the end of the period
extrarows1$ageend <- as.numeric(difftime(as.Date("1839-12-31"),extrarows1$bdate,units="days") / 365.25)
extrarows2$ageend <- as.numeric(difftime(as.Date("1869-12-31"),extrarows2$bdate,units="days") / 365.25)
extrarows3$ageend <- as.numeric(difftime(as.Date("1889-12-31"),extrarows3$bdate,units="days") / 365.25)
extrarows4$ageend <- as.numeric(difftime(as.Date("1912-12-31"),extrarows4$bdate,units="days") / 365.25)

extrarows1$end <- extrarows1$ageend
extrarows2$end <- extrarows2$ageend
extrarows3$end <- extrarows3$ageend
extrarows4$end <- extrarows4$ageend

#unless the last observation age is actually earlier, indicated by age at last observation
extrarows1$llage <- as.numeric(difftime(extrarows1$lldate, extrarows1$bdate,units="days") / 365.25)
extrarows2$llage <- as.numeric(difftime(extrarows2$lldate, extrarows2$bdate,units="days") / 365.25)
extrarows3$llage <- as.numeric(difftime(extrarows3$lldate, extrarows3$bdate,units="days") / 365.25)
extrarows4$llage <- as.numeric(difftime(extrarows4$lldate, extrarows4$bdate,units="days") / 365.25)

extrarows1$end[extrarows1$ageend > extrarows1$agell & extrarows1$lldate < as.Date("1840-01-01")] <- extrarows1$llage[extrarows1$ageend > extrarows1$llage & extrarows1$lldate < as.Date("1840-01-01")] 
extrarows2$end[extrarows2$ageend > extrarows2$agell & extrarows2$lldate < as.Date("1870-01-01")] <- extrarows2$llage[extrarows2$ageend > extrarows2$llage & extrarows2$lldate < as.Date("1870-01-01")] 
extrarows3$end[extrarows3$ageend > extrarows3$agell & extrarows3$lldate < as.Date("1890-01-01")] <- extrarows3$llage[extrarows3$ageend > extrarows3$llage & extrarows3$lldate < as.Date("1890-01-01")] 
extrarows4$end[extrarows4$ageend > extrarows4$agell & extrarows4$lldate < as.Date("1913-01-01")] <- extrarows4$llage[extrarows4$ageend > extrarows4$llage & extrarows4$lldate < as.Date("1913-01-01")] 

#if end age is over 45 within this period, censor at age 45
extrarows1$end[extrarows1$ageend > 45] <- 45
extrarows2$end[extrarows2$ageend > 45] <- 45
extrarows3$end[extrarows3$ageend > 45] <- 45
extrarows4$end[extrarows4$ageend > 45] <- 45

#if it's a death before age 45 and before end of this interval, update to age at death
extrarows1$end[extrarows1$aged < 45 & extrarows1$dyr < 1840 & extrarows1$aged > 20 & extrarows1$dyr > 1811] <- extrarows1$aged[extrarows1$aged < 45 & extrarows1$dyr < 1840 & extrarows1$aged > 20 & extrarows1$dyr > 1811]
extrarows2$end[extrarows2$aged < 45 & extrarows2$dyr < 1870 & extrarows2$aged > 20 & extrarows2$dyr > 1839] <- extrarows2$aged[extrarows2$aged < 45 & extrarows2$dyr < 1870 & extrarows2$aged > 20 & extrarows2$dyr > 1839]
extrarows3$end[extrarows3$aged < 45 & extrarows3$dyr < 1890 & extrarows3$aged > 20 & extrarows3$dyr > 1869] <- extrarows3$aged[extrarows3$aged < 45 & extrarows3$dyr < 1890 & extrarows3$aged > 20 & extrarows3$dyr > 1869]
extrarows4$end[extrarows4$aged < 45 & extrarows4$dyr < 1913 & extrarows4$aged > 20 & extrarows4$dyr > 1889] <- extrarows4$aged[extrarows4$aged < 45 & extrarows4$dyr < 1913 & extrarows4$aged > 20 & extrarows4$dyr > 1889]

#remove those who're not in the right age range 
extrarows1 <- subset(extrarows1, extrarows1$end > 20 | extrarows1$end == 20)
extrarows2 <- subset(extrarows2, extrarows2$end > 20 | extrarows2$end == 20)
extrarows3 <- subset(extrarows3, extrarows3$end > 20 | extrarows3$end == 20)
extrarows4 <- subset(extrarows4, extrarows4$end > 20 | extrarows4$end == 20)

extrarows1 <- subset(extrarows1, extrarows1$startage < 45)
extrarows2 <- subset(extrarows2, extrarows2$startage < 45)
extrarows3 <- subset(extrarows3, extrarows3$startage < 45)
extrarows4 <- subset(extrarows4, extrarows4$startage < 45)

#event variable
extrarows1$event <- 0
extrarows2$event <- 0
extrarows3$event <- 0
extrarows4$event <- 0
extrarows1$event[extrarows1$aged==extrarows1$end] <- 1
extrarows2$event[extrarows2$aged==extrarows2$end] <- 1
extrarows3$event[extrarows3$aged==extrarows3$end] <- 1
extrarows4$event[extrarows4$aged==extrarows4$end] <- 1
table(extrarows1$event) #2779 deaths in the first period
table(extrarows2$event) #5443 deaths in the 2nd period
table(extrarows3$event) #3020 deaths in the 3rd period
table(extrarows4$event) #2166 deaths in the 4th period

#add here a row where all the data rows are copied below each other into a dataframe called anl2
anl22 <- rbind(extrarows1, extrarows2, extrarows3, extrarows4)
table(anl22$event)
length(unique(anl22$personid))
anl2 <- anl22

#make a vector of all the mother ids
vecmothers2 <- unique(anl2$personid)

#a bit more data cleaning in the child sample
anl <- subset(anl, anl$mother %in% vecmothers2)

#check the total number, reported in appendix figure A1
#mothers
length(unique(anl$mother))
#children
length(unique(anl$personid))
#events
table(anl$event)
table(anl$bcohort, anl$event)

#add the highest observed paternal social class of women's children
anl2$his5 <- NULL
ses <- persons[,c("mother", "his5")]
ses <- ses[with(ses, order(ses$mother, ses$his5)),]
ses <- ses[-which(duplicated(ses$mother)),]
anl2 <- merge(anl2, ses, by.x="personid", by.y="mother", all.x=T)
#set reference category
anl2$his5 <- factor(anl2$his5, levels=c("5. Unskilled workers", "1. Elite", "2. Lower middle class", "3. Skilled workers", "4. Farmers and fishermen", "0. Missing"))
str(anl2)

#descriptive table
#mother sample
#check the total number mothers
length(unique(anl2$personid))
#events
table(anl2$event)
table(anl2$bcohort, anl2$event)
table(anl2$bcohort)
test <- subset(anl2, anl2$bcohort==1)
length(unique(test$personid))
test <- subset(anl2, anl2$bcohort==2)
length(unique(test$personid))
test <- subset(anl2, anl2$bcohort==3)
length(unique(test$personid))
test <- subset(anl2, anl2$bcohort==4)
length(unique(test$personid))

length(unique(anl2$personid)) == length(unique(anl$mother)) #true

#ses descriptives of the mothers
anl2$exposure <- anl2$end - anl2$start
sum(anl2$exposure)
table(anl2$his5)
#descriptive table, mother sample, column 1
round(((sum(anl2$exposure[anl2$his5=="1. Elite"]) / sum(anl2$exposure)) * 100), 1)
round(((sum(anl2$exposure[anl2$his5=="2. Lower middle class"]) / sum(anl2$exposure)) * 100), 1)
round(((sum(anl2$exposure[anl2$his5=="3. Skilled workers"]) / sum(anl2$exposure)) * 100), 1)
round(((sum(anl2$exposure[anl2$his5=="4. Farmers and fishermen"]) / sum(anl2$exposure)) * 100), 1)
round((sum(anl2$exposure[anl2$his5=="5. Unskilled workers"]) / sum(anl2$exposure) * 100), 1)
round(((sum(anl2$exposure[anl2$his5=="0. Missing"]) / sum(anl2$exposure) * 100)), 1)
#subset by cohort
#redo for each birth cohort
cohort <- function(){
a <- round(((sum(cohort1$exposure[cohort1$his5=="1. Elite"]) / sum(cohort1$exposure)) * 100), 1)
b <- round(((sum(cohort1$exposure[cohort1$his5=="2. Lower middle class"]) / sum(cohort1$exposure)) * 100), 1)
d <- round(((sum(cohort1$exposure[cohort1$his5=="3. Skilled workers"]) / sum(cohort1$exposure)) * 100), 1)
e <- round(((sum(cohort1$exposure[cohort1$his5=="4. Farmers and fishermen"]) / sum(cohort1$exposure)) * 100), 1)
f <- round((sum(cohort1$exposure[cohort1$his5=="5. Unskilled workers"]) / sum(cohort1$exposure) * 100), 1)
g <- round(((sum(cohort1$exposure[cohort1$his5=="0. Missing"]) / sum(cohort1$exposure) * 100)), 1)
h <- sum(cohort1$exposure)
print(c(a,b,d,e,f,g,h))
}
cohort1 <- subset(anl2, anl2$bcohort==1)
cohort()
cohort1 <- subset(anl2, anl2$bcohort==2)
cohort()
cohort1 <- subset(anl2, anl2$bcohort==3)
cohort()
cohort1 <- subset(anl2, anl2$bcohort==4)
cohort()

#descriptives child sample
table(anl$event)
length(unique(anl$mother))
table(anl$bcohort, anl$event)
table(anl$bcohort)
test <- subset(anl, anl$bcohort==1)
length(unique(test$personid))
length(unique(test$mother))
test <- subset(anl, anl$bcohort==2)
length(unique(test$personid))
length(unique(test$mother))
test <- subset(anl, anl$bcohort==3)
length(unique(test$personid))
length(unique(test$mother))
test <- subset(anl, anl$bcohort==4)
length(unique(test$personid))
length(unique(test$mother))

round(table(anl$his5) / length(unique(anl$personid)) * 100,1)
#his5 by birth cohort
table(anl$his5, anl$bcohort)
table(anl$his5, anl$bcohort)[,1]
round(table(anl$his5, anl$bcohort)[,1] / as.numeric(table(anl$bcohort)[1]) * 100,1)
table(anl$his5, anl$bcohort)[,2]
round(table(anl$his5, anl$bcohort)[,2] / as.numeric(table(anl$bcohort)[2]) * 100,1)
table(anl$his5, anl$bcohort)[,3]
round(table(anl$his5, anl$bcohort)[,3] / as.numeric(table(anl$bcohort)[3]) * 100,1)
table(anl$his5, anl$bcohort)[,4]
round(table(anl$his5, anl$bcohort)[,4] / as.numeric(table(anl$bcohort)[4]) * 100,1)


#check SES by bcohort, reported in descriptive table
table(anl$his5)
round(table(anl$his5) / length(unique(anl$personid)) * 100,1)
round(table(anl$his5, anl$bcohort) / as.numeric(table(anl$bcohort)) * 100,1)
table(anl$his5, anl$event)
table(anl$bcohort)
table(anl$bcohort, anl$event)
table(anl$event)
period11 <- subset(anl, anl$bcohort==1)

#maternal mortality rate by # of child births
#should it be live child births?

mmr <- anl
str(mmr)
mmr$var1 <- 1
table(mmr$event)

#Construct the number of births in each year
Yearlybirths <- subset(mmr, aged!="0")
agg1 <- aggregate(var1 ~ byr, sum, data = Yearlybirths)
colnames(agg1)[c(1,2)] <- c("Year", "births")

#Calculate maternal deaths per year
mmr$eventyear[mmr$event==1] <- mmr$byr[mmr$event==1]
motherdeath <- subset(mmr, !is.na(mmr$eventyear))
agg2 <- aggregate(var1 ~ eventyear, sum, data = motherdeath)
colnames(agg2)[1] <- "Year"
colnames(agg2)[2] <- "motherdeaths"

#Merge all dataframes and select relevant period
agg5 <- Reduce(function(x, y) merge(x, y, all=TRUE), list(agg1, agg2))
agg5 <- subset(agg5, Year >= 1812 & Year < 1914)
agg5$mmr <- (agg5$motherdeaths / agg5$births) * 100000
mean(agg5$mmr)
write.csv(agg5, file="matmortlinks.csv")

matmort1 <- subset(matmort, matmort$where=="nederland")
matmort2 <- subset(matmort, matmort$where=="zeeland")

#draw the graph 
png(file="fig1_mmrate.png", units="in", width=10, height=6, res=300)
#maternal mortality
matmort99 <- plot(agg5$Year, agg5$mmr, xlab="Year", ylab="Maternal mortality per 100.000 live births", type="l", lty=1,cex.lab=1.1, lwd=2, ylim=c(0,2000), xlim=c(1813, 1915), xaxs = "i", yaxs = "i", bty="n", col="pink")
#national
matmort99 <- lines(matmort1$year, matmort1$mmr, pch=23, lty=1, lwd=2,col="grey") #national
matmort99 <- lines(matmort2$year, matmort2$mmr, pch=23, lty=1, lwd=2,col="orange") #regional

axis(side=1, at = c(1812, seq(1820, 1910, by=5), 1913))
legend(x=1815, y=400, legend = c("From sample", "National", "Zeeland"), lwd=2.3, bty="n", col=c("pink", "grey", "orange")) 
matmort1
dev.off()
getwd()



############child sample
#####ses differences by parental social class

period1 <- anl

table2 <- function(i){
  j <- 1
  period1 <- subset(anl, period1$bcohort==i)
  table(period1$event)
  length(unique(period1$personid))
  length(unique(period1$mother))
  table(period1$his5)

  #ses differences
  a1 <- coxph(Surv(start, end, event) ~ his5, data=period1)
  print(a1)
  
  #estimate
  ses1 <- round(exp(a1$coefficients)[1], 3)
  ses2 <- round(exp(a1$coefficients)[2], 3)
  ses3 <- round(exp(a1$coefficients)[3], 3)
  ses4 <- round(exp(a1$coefficients)[4], 3)
  ses5 <- round(exp(a1$coefficients)[5], 3)

  #confidence intervals, lower
  print(exp(confint(a1)))
  ses1b <- round(exp(confint(a1))[1,1], 3)
  ses2b <- round(exp(confint(a1))[2,1], 3)
  ses3b <- round(exp(confint(a1))[3,1], 3)
  ses4b <- round(exp(confint(a1))[4,1], 3)
  ses5b <- round(exp(confint(a1))[5,1], 3)
 
  #confidence intervals, upper
  print(exp(confint(a1)))
  ses1c <- round(exp(confint(a1))[1,2], 3)
  ses2c <- round(exp(confint(a1))[2,2], 3)
  ses3c <- round(exp(confint(a1))[3,2], 3)
  ses4c <- round(exp(confint(a1))[4,2], 3)
  ses5c <- round(exp(confint(a1))[5,2], 3)
 
  row <- c(i,j,length(unique(period1$personid)),a1$nevent,ses1, ses1b, ses1c, ses2, ses2b, ses2c, ses3, ses3b, ses3c, ses4, ses4b, ses4c, ses5, ses5b, ses5c)
  row1 <<- rbind(row1, row)
}

row1 <- c("year1", "year2", "N individuals", "N Deaths", "SES elite", "SES elite CI low", "SES elite CI high", 
          "SES lower middle class", "SES lower middle class CI low", "SES lower middle class CI high", 
          "SES skilled", "SES skilled CI low", "SES skilled CI high", 
          "farmers", "farmers CI low","farmers CI high",
          "missing", "missing CI low", "missing CI high")
table2(1)
table2(2)
table2(3)
table2(4)
write.csv(row1, "figure2.csv") 

#figure
table2year <- read.csv("figure2.csv")
table2year <- table2year[-1,-c(1:3)]
colnames(table2year) <- c("N", "Deaths", "SES1", "SES1CI1", "SES1CI2", "SES2", "SES2CI1", "SES2CI2", 
                          "SES3", "SES3CI1", "SES3CI2", "SES4", "SES4CI1","SES4CI2",
                          "SES5", "SES5CI1", "SES5CI2")
table2year <- as.data.frame(lapply(table2year, as.numeric))
table2year$Period <- c("1812-1839", "1840-1869", "1870-1889", "1890-1913")


#Plot 1: higher managers
label <- table2year$Period
mean <- table2year$SES1
lower <- table2year$SES1CI1
upper <- table2year$SES1CI2
df <- data.frame(label, mean, lower, upper)
df

#the plot
df$label <- factor(df$label, levels=rev(df$label)) #reverse labels
fp2 <- ggplot(data=df, aes(x=label, y=mean, ymin=lower, ymax=upper)) +
  geom_pointrange() + 
  geom_hline(yintercept=1, lty=2, col="darkred") +  #add a dotted line at x=1 after flip
  coord_flip() +  # flip coordinates (puts labels on y axis)
  xlab("") + ylab("") +
  theme_bw()

#adjust the looks
fp2 <- fp2 + ylim(0.5,2)
fp2 <- fp2 + theme_bw()
fp2 <- fp2 + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
fp2 <- fp2 + theme(text = element_text(size=12, family = 'sans'))
fp2 <- fp2 + ggtitle(label="Elite")
plot1 <- fp2; plot1

#plot2
label <- table2year$Period
mean <- table2year$SES2
lower <- table2year$SES2CI1
upper <- table2year$SES2CI2
df <- data.frame(label, mean, lower, upper)

#the plot
df$label <- factor(df$label, levels=rev(df$label)) #reverse labels
fp2 <- ggplot(data=df, aes(x=label, y=mean, ymin=lower, ymax=upper)) +
  geom_pointrange() + 
  geom_hline(yintercept=1, lty=2, col="darkred") +  #add a dotted line at x=1 after flip
  coord_flip() +  # flip coordinates (puts labels on y axis)
  xlab("") + ylab("") +
  theme_bw()

#adjust the looks
fp2 <- fp2 + ylim(0.5,2)
fp2 <- fp2 + theme_bw()
fp2 <- fp2 + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
fp2 <- fp2 + theme(text = element_text(size=12, family = 'sans'))
fp2 <- fp2 + ggtitle(label="Lower middle class")
plot2 <- fp2; plot2

#plot3
label <- table2year$Period
mean <- table2year$SES3
lower <- table2year$SES3CI1
upper <- table2year$SES3CI2
df <- data.frame(label, mean, lower, upper)

#the plot
df$label <- factor(df$label, levels=rev(df$label)) #reverse labels
fp2 <- ggplot(data=df, aes(x=label, y=mean, ymin=lower, ymax=upper)) +
  geom_pointrange() + 
  geom_hline(yintercept=1, lty=2, col="darkred") +  #add a dotted line at x=1 after flip
  coord_flip() +  # flip coordinates (puts labels on y axis)
  xlab("") + ylab("") +
  theme_bw()

#adjust the looks
fp2 <- fp2 + ylim(0.5,2)
fp2 <- fp2 + theme_bw()
fp2 <- fp2 + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
fp2 <- fp2 + theme(text = element_text(size=12, family = 'sans'))
fp2 <- fp2 + ggtitle(label="Skilled workers")
plot3 <- fp2; plot3


#plot4
label <- table2year$Period
mean <- table2year$SES4
lower <- table2year$SES4CI1
upper <- table2year$SES4CI2
df <- data.frame(label, mean, lower, upper)

#the plot
df$label <- factor(df$label, levels=rev(df$label)) #reverse labels
fp2 <- ggplot(data=df, aes(x=label, y=mean, ymin=lower, ymax=upper)) +
  geom_pointrange() + 
  geom_hline(yintercept=1, lty=2, col="darkred") +  #add a dotted line at x=1 after flip
  coord_flip() +  # flip coordinates (puts labels on y axis)
  xlab("") + ylab("") +
  theme_bw()

#adjust the looks
fp2 <- fp2 + ylim(0.5,2)
fp2 <- fp2 + theme_bw()
fp2 <- fp2 + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
fp2 <- fp2 + theme(text = element_text(size=12, family = 'sans'))
fp2 <- fp2 + ggtitle(label="Farmers")
plot4 <- fp2; plot4

plot9 <- ggarrange(plot1, plot2, plot3, plot4, ncol=4, nrow=1)
plot9
ggsave(plot=plot9, filename="fig2_results_his5.pdf", width=16, height=4, units="in", dpi=300) 
#this figure is figure 2



################general mortality of women by social class

period1 <- anl2

table2 <- function(i){
  j <- 1
  period1 <- subset(period1, period1$bcohort==i)
  str(period1)
  table(period1$event)
  length(unique(period1$personid))
  table(period1$his5)
  
  #ses differences
  a1 <- coxph(Surv(startage, end, event) ~ his5, data=period1)
  print(a1)
  
  #estimate
  ses1 <- round(exp(a1$coefficients)[1], 3)
  ses2 <- round(exp(a1$coefficients)[2], 3)
  ses3 <- round(exp(a1$coefficients)[3], 3)
  ses4 <- round(exp(a1$coefficients)[4], 3)
  ses5 <- round(exp(a1$coefficients)[5], 3)
  
  #confidence intervals, lower
  print(exp(confint(a1)))
  ses1b <- round(exp(confint(a1))[1,1], 3)
  ses2b <- round(exp(confint(a1))[2,1], 3)
  ses3b <- round(exp(confint(a1))[3,1], 3)
  ses4b <- round(exp(confint(a1))[4,1], 3)
  ses5b <- round(exp(confint(a1))[5,1], 3)
  
  #confidence intervals, upper
  print(exp(confint(a1)))
  ses1c <- round(exp(confint(a1))[1,2], 3)
  ses2c <- round(exp(confint(a1))[2,2], 3)
  ses3c <- round(exp(confint(a1))[3,2], 3)
  ses4c <- round(exp(confint(a1))[4,2], 3)
  ses5c <- round(exp(confint(a1))[5,2], 3)
  
  row <- c(i,j,length(unique(period1$personid)),a1$nevent,ses1, ses1b, ses1c, ses2, ses2b, ses2c, ses3, ses3b, ses3c, ses4, ses4b, ses4c, ses5, ses5b, ses5c)
  row1 <<- rbind(row1, row)
}

row1 <- c("year1", "year2", "N individuals", "N Deaths", "SES elite", "SES elite CI low", "SES elite CI high", 
          "SES lower middle class", "SES lower middle class CI low", "SES lower middle class CI high", 
          "SES skilled", "SES skilled CI low", "SES skilled CI high", 
          "farmers", "farmers CI low","farmers CI high",
          "missing", "missing CI low", "missing CI high")
table2(1)
table2(2)
table2(3)
table2(4)
write.csv(row1, "figure3.csv") 

#figure
table2year <- read.csv("figure3.csv")
table2year <- table2year[-1,-c(1:3)]
colnames(table2year) <- c("N", "Deaths", "SES1", "SES1CI1", "SES1CI2", "SES2", "SES2CI1", "SES2CI2", 
                          "SES3", "SES3CI1", "SES3CI2", "SES4", "SES4CI1","SES4CI2",
                          "SES5", "SES5CI1", "SES5CI2")
table2year <- as.data.frame(lapply(table2year, as.numeric))
table2year$Period <- c("1812-1839", "1840-1869", "1870-1889", "1890-1913")


#Plot 1: higher managers
label <- table2year$Period
mean <- table2year$SES1
lower <- table2year$SES1CI1
upper <- table2year$SES1CI2
df <- data.frame(label, mean, lower, upper)
df

#the plot
df$label <- factor(df$label, levels=rev(df$label)) #reverse labels
fp2 <- ggplot(data=df, aes(x=label, y=mean, ymin=lower, ymax=upper)) +
  geom_pointrange() + 
  geom_hline(yintercept=1, lty=2, col="darkred") +  #add a dotted line at x=1 after flip
  coord_flip() +  # flip coordinates (puts labels on y axis)
  xlab("") + ylab("") +
  theme_bw()

#adjust the looks
fp2 <- fp2 + ylim(0.5,2)
fp2 <- fp2 + theme_bw()
fp2 <- fp2 + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
fp2 <- fp2 + theme(text = element_text(size=12, family = 'sans'))
fp2 <- fp2 + ggtitle(label="Elite")
plot1 <- fp2; plot1

#plot2
label <- table2year$Period
mean <- table2year$SES2
lower <- table2year$SES2CI1
upper <- table2year$SES2CI2
df <- data.frame(label, mean, lower, upper)

#the plot
df$label <- factor(df$label, levels=rev(df$label)) #reverse labels
fp2 <- ggplot(data=df, aes(x=label, y=mean, ymin=lower, ymax=upper)) +
  geom_pointrange() + 
  geom_hline(yintercept=1, lty=2, col="darkred") +  #add a dotted line at x=1 after flip
  coord_flip() +  # flip coordinates (puts labels on y axis)
  xlab("") + ylab("") +
  theme_bw()

#adjust the looks
fp2 <- fp2 + ylim(0.5,2)
fp2 <- fp2 + theme_bw()
fp2 <- fp2 + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
fp2 <- fp2 + theme(text = element_text(size=12, family = 'sans'))
fp2 <- fp2 + ggtitle(label="Lower middle class")
plot2 <- fp2; plot2

#plot3
label <- table2year$Period
mean <- table2year$SES3
lower <- table2year$SES3CI1
upper <- table2year$SES3CI2
df <- data.frame(label, mean, lower, upper)

#the plot
df$label <- factor(df$label, levels=rev(df$label)) #reverse labels
fp2 <- ggplot(data=df, aes(x=label, y=mean, ymin=lower, ymax=upper)) +
  geom_pointrange() + 
  geom_hline(yintercept=1, lty=2, col="darkred") +  #add a dotted line at x=1 after flip
  coord_flip() +  # flip coordinates (puts labels on y axis)
  xlab("") + ylab("") +
  theme_bw()

#adjust the looks
fp2 <- fp2 + ylim(0.5,2)
fp2 <- fp2 + theme_bw()
fp2 <- fp2 + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
fp2 <- fp2 + theme(text = element_text(size=12, family = 'sans'))
fp2 <- fp2 + ggtitle(label="Skilled workers")
plot3 <- fp2; plot3


#plot4
label <- table2year$Period
mean <- table2year$SES4
lower <- table2year$SES4CI1
upper <- table2year$SES4CI2
df <- data.frame(label, mean, lower, upper)

#the plot
df$label <- factor(df$label, levels=rev(df$label)) #reverse labels
fp2 <- ggplot(data=df, aes(x=label, y=mean, ymin=lower, ymax=upper)) +
  geom_pointrange() + 
  geom_hline(yintercept=1, lty=2, col="darkred") +  #add a dotted line at x=1 after flip
  coord_flip() +  # flip coordinates (puts labels on y axis)
  xlab("") + ylab("") +
  theme_bw()

#adjust the looks
fp2 <- fp2 + ylim(0.5,2)
fp2 <- fp2 + theme_bw()
fp2 <- fp2 + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
fp2 <- fp2 + theme(text = element_text(size=12, family = 'sans'))
fp2 <- fp2 + ggtitle(label="Farmers")
plot4 <- fp2; plot4

plot9 <- ggarrange(plot1, plot2, plot3, plot4, ncol=4, nrow=1)
plot9
ggsave(plot=plot9, filename="fig3_results_his5_allmort.pdf", width=16, height=4, units="in", dpi=300) 


#same but without maternal mortality included

period1 <- anl2
#redefine event: =0 if in the first 42 days
str(period1)
maternaldeaths <- anl[,c("mother", "matmort", "ddatem")]
maternaldeaths <- subset(maternaldeaths, matmort==1)
colnames(maternaldeaths) <- c("personid", "matmort2", "ddatem")
period1 <- merge(period1, maternaldeaths, by="personid", all.x=T)
period1$event2 <- period1$event #all 13,408 events
table(period1$event2)
period1$event2[period1$matmort2==1] <- 0
table(period1$event2)

table2 <- function(i){
  j <- 1
  period1 <- subset(period1, period1$bcohort==i)
  str(period1)
  table(period1$event2)
  length(unique(period1$personid))
  table(period1$his5)
  
  #ses differences
  a1 <- coxph(Surv(startage, end, event2) ~ his5, data=period1)
  print(a1)
  
  #estimate
  ses1 <- round(exp(a1$coefficients)[1], 3)
  ses2 <- round(exp(a1$coefficients)[2], 3)
  ses3 <- round(exp(a1$coefficients)[3], 3)
  ses4 <- round(exp(a1$coefficients)[4], 3)
  ses5 <- round(exp(a1$coefficients)[5], 3)
  
  #confidence intervals, lower
  print(exp(confint(a1)))
  ses1b <- round(exp(confint(a1))[1,1], 3)
  ses2b <- round(exp(confint(a1))[2,1], 3)
  ses3b <- round(exp(confint(a1))[3,1], 3)
  ses4b <- round(exp(confint(a1))[4,1], 3)
  ses5b <- round(exp(confint(a1))[5,1], 3)
  
  #confidence intervals, upper
  print(exp(confint(a1)))
  ses1c <- round(exp(confint(a1))[1,2], 3)
  ses2c <- round(exp(confint(a1))[2,2], 3)
  ses3c <- round(exp(confint(a1))[3,2], 3)
  ses4c <- round(exp(confint(a1))[4,2], 3)
  ses5c <- round(exp(confint(a1))[5,2], 3)
  
  row <- c(i,j,length(unique(period1$personid)),a1$nevent,ses1, ses1b, ses1c, ses2, ses2b, ses2c, ses3, ses3b, ses3c, ses4, ses4b, ses4c, ses5, ses5b, ses5c)
  row1 <<- rbind(row1, row)
}

row1 <- c("year1", "year2", "N individuals", "N Deaths", "SES elite", "SES elite CI low", "SES elite CI high", 
          "SES lower middle class", "SES lower middle class CI low", "SES lower middle class CI high", 
          "SES skilled", "SES skilled CI low", "SES skilled CI high", 
          "farmers", "farmers CI low","farmers CI high",
          "missing", "missing CI low", "missing CI high")
table2(1)
table2(2)
table2(3)
table2(4)
write.csv(row1, "figure4.csv") 

#figure
table2year <- read.csv("figure4.csv")
table2year <- table2year[-1,-c(1:3)]
colnames(table2year) <- c("N", "Deaths", "SES1", "SES1CI1", "SES1CI2", "SES2", "SES2CI1", "SES2CI2", 
                          "SES3", "SES3CI1", "SES3CI2", "SES4", "SES4CI1","SES4CI2",
                          "SES5", "SES5CI1", "SES5CI2")
table2year <- as.data.frame(lapply(table2year, as.numeric))
table2year$Period <- c("1812-1839", "1840-1869", "1870-1889", "1890-1913")


#Plot 1: higher managers
label <- table2year$Period
mean <- table2year$SES1
lower <- table2year$SES1CI1
upper <- table2year$SES1CI2
df <- data.frame(label, mean, lower, upper)
df

#the plot
df$label <- factor(df$label, levels=rev(df$label)) #reverse labels
fp2 <- ggplot(data=df, aes(x=label, y=mean, ymin=lower, ymax=upper)) +
  geom_pointrange() + 
  geom_hline(yintercept=1, lty=2, col="darkred") +  #add a dotted line at x=1 after flip
  coord_flip() +  # flip coordinates (puts labels on y axis)
  xlab("") + ylab("") +
  theme_bw()

#adjust the looks
fp2 <- fp2 + ylim(0.5,2)
fp2 <- fp2 + theme_bw()
fp2 <- fp2 + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
fp2 <- fp2 + theme(text = element_text(size=12, family = 'sans'))
fp2 <- fp2 + ggtitle(label="Elite")
plot1 <- fp2; plot1

#plot2
label <- table2year$Period
mean <- table2year$SES2
lower <- table2year$SES2CI1
upper <- table2year$SES2CI2
df <- data.frame(label, mean, lower, upper)

#the plot
df$label <- factor(df$label, levels=rev(df$label)) #reverse labels
fp2 <- ggplot(data=df, aes(x=label, y=mean, ymin=lower, ymax=upper)) +
  geom_pointrange() + 
  geom_hline(yintercept=1, lty=2, col="darkred") +  #add a dotted line at x=1 after flip
  coord_flip() +  # flip coordinates (puts labels on y axis)
  xlab("") + ylab("") +
  theme_bw()

#adjust the looks
fp2 <- fp2 + ylim(0.5,2)
fp2 <- fp2 + theme_bw()
fp2 <- fp2 + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
fp2 <- fp2 + theme(text = element_text(size=12, family = 'sans'))
fp2 <- fp2 + ggtitle(label="Lower middle class")
plot2 <- fp2; plot2

#plot3
label <- table2year$Period
mean <- table2year$SES3
lower <- table2year$SES3CI1
upper <- table2year$SES3CI2
df <- data.frame(label, mean, lower, upper)

#the plot
df$label <- factor(df$label, levels=rev(df$label)) #reverse labels
fp2 <- ggplot(data=df, aes(x=label, y=mean, ymin=lower, ymax=upper)) +
  geom_pointrange() + 
  geom_hline(yintercept=1, lty=2, col="darkred") +  #add a dotted line at x=1 after flip
  coord_flip() +  # flip coordinates (puts labels on y axis)
  xlab("") + ylab("") +
  theme_bw()

#adjust the looks
fp2 <- fp2 + ylim(0.5,2)
fp2 <- fp2 + theme_bw()
fp2 <- fp2 + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
fp2 <- fp2 + theme(text = element_text(size=12, family = 'sans'))
fp2 <- fp2 + ggtitle(label="Skilled workers")
plot3 <- fp2; plot3


#plot4
label <- table2year$Period
mean <- table2year$SES4
lower <- table2year$SES4CI1
upper <- table2year$SES4CI2
df <- data.frame(label, mean, lower, upper)

#the plot
df$label <- factor(df$label, levels=rev(df$label)) #reverse labels
fp2 <- ggplot(data=df, aes(x=label, y=mean, ymin=lower, ymax=upper)) +
  geom_pointrange() + 
  geom_hline(yintercept=1, lty=2, col="darkred") +  #add a dotted line at x=1 after flip
  coord_flip() +  # flip coordinates (puts labels on y axis)
  xlab("") + ylab("") +
  theme_bw()

#adjust the looks
fp2 <- fp2 + ylim(0.5,2)
fp2 <- fp2 + theme_bw()
fp2 <- fp2 + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
fp2 <- fp2 + theme(text = element_text(size=12, family = 'sans'))
fp2 <- fp2 + ggtitle(label="Farmers")
plot4 <- fp2; plot4

plot9 <- ggarrange(plot1, plot2, plot3, plot4, ncol=4, nrow=1)
plot9
ggsave(plot=plot9, filename="fig4_results_his5_allmortbutnomatmort.pdf", width=16, height=4, units="in", dpi=300) 






##results for only women for whom a death date or other follow up is known in Zeeland, for social class models
str(anl22)
anl23 <- subset(anl2, anl2$llage > 45 | anl2$event==1)


period1 <- anl23

table2 <- function(i){
  j <- 1
  period1 <- subset(period1, period1$bcohort==i)
  str(period1)
  table(period1$event)
  length(unique(period1$personid))
  table(period1$his5)
  
  #ses differences
  a1 <- coxph(Surv(startage, end, event) ~ his5, data=period1)
  print(a1)
  
  #estimate
  ses1 <- round(exp(a1$coefficients)[1], 3)
  ses2 <- round(exp(a1$coefficients)[2], 3)
  ses3 <- round(exp(a1$coefficients)[3], 3)
  ses4 <- round(exp(a1$coefficients)[4], 3)
  ses5 <- round(exp(a1$coefficients)[5], 3)
  
  #confidence intervals, lower
  print(exp(confint(a1)))
  ses1b <- round(exp(confint(a1))[1,1], 3)
  ses2b <- round(exp(confint(a1))[2,1], 3)
  ses3b <- round(exp(confint(a1))[3,1], 3)
  ses4b <- round(exp(confint(a1))[4,1], 3)
  ses5b <- round(exp(confint(a1))[5,1], 3)
  
  #confidence intervals, upper
  print(exp(confint(a1)))
  ses1c <- round(exp(confint(a1))[1,2], 3)
  ses2c <- round(exp(confint(a1))[2,2], 3)
  ses3c <- round(exp(confint(a1))[3,2], 3)
  ses4c <- round(exp(confint(a1))[4,2], 3)
  ses5c <- round(exp(confint(a1))[5,2], 3)
  
  row <- c(i,j,length(unique(period1$personid)),a1$nevent,ses1, ses1b, ses1c, ses2, ses2b, ses2c, ses3, ses3b, ses3c, ses4, ses4b, ses4c, ses5, ses5b, ses5c)
  row1 <<- rbind(row1, row)
}

row1 <- c("year1", "year2", "N individuals", "N Deaths", "SES elite", "SES elite CI low", "SES elite CI high", 
          "SES lower middle class", "SES lower middle class CI low", "SES lower middle class CI high", 
          "SES skilled", "SES skilled CI low", "SES skilled CI high", 
          "farmers", "farmers CI low","farmers CI high",
          "missing", "missing CI low", "missing CI high")
table2(1)
table2(2)
table2(3)
table2(4)
write.csv(row1, "figureA2.csv") 

#figure
table2year <- read.csv("figureA2.csv")
table2year <- table2year[-1,-c(1:3)]
colnames(table2year) <- c("N", "Deaths", "SES1", "SES1CI1", "SES1CI2", "SES2", "SES2CI1", "SES2CI2", 
                          "SES3", "SES3CI1", "SES3CI2", "SES4", "SES4CI1","SES4CI2",
                          "SES5", "SES5CI1", "SES5CI2")
table2year <- as.data.frame(lapply(table2year, as.numeric))
table2year$Period <- c("1812-1839", "1840-1869", "1870-1889", "1890-1913")


#Plot 1: higher managers
label <- table2year$Period
mean <- table2year$SES1
lower <- table2year$SES1CI1
upper <- table2year$SES1CI2
df <- data.frame(label, mean, lower, upper)
df

#the plot
df$label <- factor(df$label, levels=rev(df$label)) #reverse labels
fp2 <- ggplot(data=df, aes(x=label, y=mean, ymin=lower, ymax=upper)) +
  geom_pointrange() + 
  geom_hline(yintercept=1, lty=2, col="darkred") +  #add a dotted line at x=1 after flip
  coord_flip() +  # flip coordinates (puts labels on y axis)
  xlab("") + ylab("") +
  theme_bw()

#adjust the looks
fp2 <- fp2 + ylim(0.5,2)
fp2 <- fp2 + theme_bw()
fp2 <- fp2 + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
fp2 <- fp2 + theme(text = element_text(size=12, family = 'sans'))
fp2 <- fp2 + ggtitle(label="Elite")
plot1 <- fp2; plot1

#plot2
label <- table2year$Period
mean <- table2year$SES2
lower <- table2year$SES2CI1
upper <- table2year$SES2CI2
df <- data.frame(label, mean, lower, upper)

#the plot
df$label <- factor(df$label, levels=rev(df$label)) #reverse labels
fp2 <- ggplot(data=df, aes(x=label, y=mean, ymin=lower, ymax=upper)) +
  geom_pointrange() + 
  geom_hline(yintercept=1, lty=2, col="darkred") +  #add a dotted line at x=1 after flip
  coord_flip() +  # flip coordinates (puts labels on y axis)
  xlab("") + ylab("") +
  theme_bw()

#adjust the looks
fp2 <- fp2 + ylim(0.5,2)
fp2 <- fp2 + theme_bw()
fp2 <- fp2 + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
fp2 <- fp2 + theme(text = element_text(size=12, family = 'sans'))
fp2 <- fp2 + ggtitle(label="Lower middle class")
plot2 <- fp2; plot2

#plot3
label <- table2year$Period
mean <- table2year$SES3
lower <- table2year$SES3CI1
upper <- table2year$SES3CI2
df <- data.frame(label, mean, lower, upper)

#the plot
df$label <- factor(df$label, levels=rev(df$label)) #reverse labels
fp2 <- ggplot(data=df, aes(x=label, y=mean, ymin=lower, ymax=upper)) +
  geom_pointrange() + 
  geom_hline(yintercept=1, lty=2, col="darkred") +  #add a dotted line at x=1 after flip
  coord_flip() +  # flip coordinates (puts labels on y axis)
  xlab("") + ylab("") +
  theme_bw()

#adjust the looks
fp2 <- fp2 + ylim(0.5,2)
fp2 <- fp2 + theme_bw()
fp2 <- fp2 + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
fp2 <- fp2 + theme(text = element_text(size=12, family = 'sans'))
fp2 <- fp2 + ggtitle(label="Skilled workers")
plot3 <- fp2; plot3


#plot4
label <- table2year$Period
mean <- table2year$SES4
lower <- table2year$SES4CI1
upper <- table2year$SES4CI2
df <- data.frame(label, mean, lower, upper)

#the plot
df$label <- factor(df$label, levels=rev(df$label)) #reverse labels
fp2 <- ggplot(data=df, aes(x=label, y=mean, ymin=lower, ymax=upper)) +
  geom_pointrange() + 
  geom_hline(yintercept=1, lty=2, col="darkred") +  #add a dotted line at x=1 after flip
  coord_flip() +  # flip coordinates (puts labels on y axis)
  xlab("") + ylab("") +
  theme_bw()

#adjust the looks
fp2 <- fp2 + ylim(0.5,2)
fp2 <- fp2 + theme_bw()
fp2 <- fp2 + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
fp2 <- fp2 + theme(text = element_text(size=12, family = 'sans'))
fp2 <- fp2 + ggtitle(label="Farmers")
plot4 <- fp2; plot4

plot9 <- ggarrange(plot1, plot2, plot3, plot4, ncol=4, nrow=1)
plot9
ggsave(plot=plot9, filename="fig_a2_results_his5_allmort_followup.pdf", width=16, height=4, units="in", dpi=300)