main="Sibshipsize of scholars by level of publications",
sub=paste("Final sample. Nobs=",nrow(profsgood),". ",pts," points, sample=",sha,sep=""),col="coral2",xaxt='n')
axis(1,at=c(1,10,20,30,40,50,60,pts),labels=c(as.integer(date[1]),date[10],date[20],date[30],date[40],date[50],date[60],date[70]))
lines(t,nbottom,lty=3,col="darkgreen",lwd=3)
legend("topright", legend=c("Above median number of publications", "Below median no. of publications"),
col=c("coral2", "darkgreen"), lty=c(1,3),lwd=3,
box.lty=0)
for (i in 1:length(sign)) {
if (p_values[i] <= 0.1) {
rect(xleft = t[i] - 0.5, xright = t[i] + 0.5, ybottom = 2.85, ytop = 3.85, col = adjustcolor("lightgreen", alpha.f = 0.3), border = NA, add = TRUE)
}
}
text(5,3.6,"evolutionary advantage")
text(5,3.57,"of the fittest")
text(65,3.3,"quality-quantity")
text(65,3.27,"tradeoff")
dev.off()
pdf(file="fertility (male) of both groups.pdf", width = 11, height = 7)
plot(t,ntop,lwd=3,type="l",xlab="Year of birth", ylab="fertility (sibshipsize male)",frame.plot=TRUE,ylim=c(2.9,3.8),
main="Sibshipsize of scholars by level of publications",
sub=paste("Final sample. Nobs=",nrow(profsgood),". ",pts," points, sample=",sha,sep=""),col="coral2",xaxt='n')
axis(1,at=c(1,10,20,30,40,50,60,pts),labels=c(as.integer(date[1]),date[10],date[20],date[30],date[40],date[50],date[60],date[70]))
lines(t,nbottom,lty=3,col="darkgreen",lwd=3)
legend("topright", legend=c("Above median number of publications", "Below median no. of publications"),
col=c("coral2", "darkgreen"), lty=c(1,3),lwd=3,
box.lty=0)
for (i in 1:length(sign)) {
if (p_values[i] <= 0.1) {
rect(xleft = t[i] - 0.5, xright = t[i] + 0.5, ybottom = 2.85, ytop = 3.85, col = adjustcolor("lightgray", alpha.f = 0.3), border = NA, add = TRUE)
}
}
text(5,3.6,"evolutionary advantage")
text(5,3.57,"of the fittest")
text(65,3.3,"quality-quantity")
text(65,3.27,"tradeoff")
dev.off()
# remove (almost) everything in the working environment.
rm(list = ls())
# to import xl files
library(readxl)
# Set working directory based on the this.path package, which reliably detects the path of the running script. and projects.
if (!require(this.path)) {
install.packages("this.path")
require(this.path)
}
setwd(dirname(this.path()))
profsgood <- as.data.frame(read_excel("siblings-extended-good.xlsx"))
#################################
profsgood <- profsgood[order(profsgood$refdate),]
seuil <- median(profsgood$titles)
seuil33 <- quantile(profsgood$titles,.33333)
seuil66 <- quantile(profsgood$titles,.66666)
profstop <- profsgood[profsgood$titles>seuil,]
profsbottom <- profsgood[profsgood$titles<=seuil,]
pts<-70 #number of points for which we compute the fertility
sha <- .3  # range of quantiles at each point
#plot of coefficient in recursive estimation
date <- (1:pts)
sign <-(1:pts)
t_test_result <- (1:pts)
p_values <- (1:pts)
t <- (1:pts)
maletop <- (1:pts)
malebottom <- (1:pts)
ns<-nrow(profsgood)
#  splitting at the median
for (i in (1:pts)){
smpl<-profsgood[as.integer(1+(i-1)*ns*(1-sha)/(pts-1)):as.integer((i-1)*ns*(1-sha)/(pts-1)+ns*sha),]
profstop <- smpl[smpl$titles>seuil,]
profsbottom <- smpl[smpl$titles<=seuil,]
maletop[i]=mean(profstop[,"MaleDescendants"],na.rm=T)
malebottom[i]=mean(profsbottom[,"MaleDescendants"],na.rm=T)
date[i]=round(mean(smpl$refdate),0)
t_test_result <- t.test(profstop[,"MaleDescendants"], profsbottom[,"MaleDescendants"])
p_values[i] <- t_test_result$p.value
}
profsgood[!is.na(profsgood$MaleDescendants) & !is.na(profsgood$NumberDescendants) & profsgood$MaleDescendants>profsgood$NumberDescendants,"PKey"]
pdf(file="background_children.pdf", width = 11, height = 7)
plot(t,maletop,lwd=3,type="l",xlab="Year of birth", ylab="fertility (male)",frame.plot=TRUE,ylim=c(1,2.5),
main="Male Children of scholars by level of publications",
sub=paste("Final sample. Nobs=",nrow(profsgood),". ",pts," points, sample=",sha,sep=""),col="coral2",xaxt='n')
axis(1,at=c(1,10,20,30,40,50,60,pts),labels=c(as.integer(date[1]),date[10],date[20],date[30],date[40],date[50],date[60],date[pts]))
lines(t,malebottom,lty=1,col="darkgreen",lwd=3)
legend("topright", legend=c("Above median number of publications", "Below median number of publications"),
col=c("coral2", "darkgreen"), lty=1,lwd=3,
box.lty=0)
#polygon(x=c(t,rev(t)),y=c(high05,rev(low05)), col =  adjustcolor("gray", alpha.f = 0.20), border = NA)
for (i in 1:length(sign)) {
if (p_values[i] > 0.1) {
rect(xleft = t[i] - 0.5, xright = t[i] + 0.5, ybottom = .9, ytop = 2.6, col = adjustcolor("lightpink", alpha.f = 0.3), border = NA, add = TRUE)
} else {
rect(xleft = t[i] - 0.5, xright = t[i] + 0.5, ybottom = .9, ytop = 2.6, col = adjustcolor("lightgreen", alpha.f = 0.3), border = NA, add = TRUE)
}
}
dev.off()
pdf(file="background_children.pdf", width = 11, height = 7)
plot(t,maletop,lwd=3,type="l",xlab="Year of birth", ylab="fertility (male)",frame.plot=TRUE,ylim=c(1,2.5),
main="Male Children of scholars by level of publications",
sub=paste("Final sample. Nobs=",nrow(profsgood),". ",pts," points, sample=",sha,sep=""),col="coral2",xaxt='n')
axis(1,at=c(1,10,20,30,40,50,60,pts),labels=c(as.integer(date[1]),date[10],date[20],date[30],date[40],date[50],date[60],date[pts]))
lines(t,malebottom,lty=3,col="darkgreen",lwd=3)
legend("topright", legend=c("Above median number of publications", "Below median number of publications"),
col=c("coral2", "darkgreen"), lty=c(1,3),lwd=3,
box.lty=0)
#polygon(x=c(t,rev(t)),y=c(high05,rev(low05)), col =  adjustcolor("gray", alpha.f = 0.20), border = NA)
for (i in 1:length(sign)) {
if (p_values[i] <= 0.1) {
rect(xleft = t[i] - 0.5, xright = t[i] + 0.5, ybottom = .9, ytop = 2.6, col = adjustcolor("lightgray", alpha.f = 0.3), border = NA, add = TRUE)
}
}
dev.off()
## remove (almost) everything in the working environment.
rm(list = ls())
# to import xl files
library(readxl)
# Set working directory based on the this.path package, which reliably detects the path of the running script. and projects.
if (!require(this.path)) {
install.packages("this.path")
require(this.path)
}
setwd(dirname(this.path()))
profsgood <- as.data.frame(read_excel("siblings-extended-good.xlsx"))
#################################
profsgood <- profsgood[order(profsgood$refdate),]
#to make it identical to regression
profsgood  <- profsgood[!is.na(profsgood$longevity),]
seuil <- median(profsgood$titles)
seuil33 <- quantile(profsgood$titles,.33333)
seuil66 <- quantile(profsgood$titles,.66666)
ns<-nrow(profsgood)
pts<-70 #number of points for which we compute the fertility
sha <- .30  # range of quantiles at each point
#plot of coefficient in recursive estimation
ntop <- (1:pts)
nbottom <- (1:pts)
date <- (1:pts)
df <- (1:pts)
t <- (1:pts)
sign <-(1:pts)
t_test_result <- (1:pts)
p_values <- (1:pts)
#dividing above/below median
for (i in (1:pts)){
smpl<-profsgood[as.integer(1+(i-1)*ns*(1-sha)/(pts-1)):as.integer((i-1)*ns*(1-sha)/(pts-1)+ns*sha),]
profstop <- smpl[smpl$titles>seuil,]
profsbottom <- smpl[smpl$titles<=seuil,]
ntop[i]=mean(profstop[,"SibshipSizeMale"])
nbottom[i]=mean(profsbottom[,"SibshipSizeMale"])
t_test_result <- t.test(profstop[,"SibshipSizeMale"], profsbottom[,"SibshipSizeMale"])
p_values[i] <- t_test_result$p.value
date[i]=round(mean(smpl$refdate),0)
}
cbind(date,ntop,nbottom)
pdf(file="fertility (male) of both groups.pdf", width = 11, height = 7)
plot(t,ntop,lwd=3,type="l",xlab="Year of birth", ylab="fertility (sibshipsize male)",frame.plot=TRUE,ylim=c(2.9,3.8),
main="Sibshipsize of scholars by level of publications",
sub=paste("Final sample. Nobs=",nrow(profsgood),". ",pts," points, sample=",sha,sep=""),col="coral2",xaxt='n')
axis(1,at=c(1,10,20,30,40,50,60,pts),labels=c(as.integer(date[1]),date[10],date[20],date[30],date[40],date[50],date[60],date[70]))
lines(t,nbottom,lty=4,col="darkgreen",lwd=3)
legend("topright", legend=c("Above median number of publications", "Below median no. of publications"),
col=c("coral2", "darkgreen"), lty=c(1,3),lwd=3,
box.lty=0)
for (i in 1:length(sign)) {
if (p_values[i] <= 0.1) {
rect(xleft = t[i] - 0.5, xright = t[i] + 0.5, ybottom = 2.85, ytop = 3.85, col = adjustcolor("lightgray", alpha.f = 0.3), border = NA, add = TRUE)
}
}
text(5,3.6,"evolutionary advantage")
text(5,3.57,"of the fittest")
text(65,3.3,"quality-quantity")
text(65,3.27,"tradeoff")
dev.off()
# remove (almost) everything in the working environment.
rm(list = ls())
# to import xl files
library(readxl)
# Set working directory based on the this.path package, which reliably detects the path of the running script. and projects.
if (!require(this.path)) {
install.packages("this.path")
require(this.path)
}
setwd(dirname(this.path()))
profsgood <- as.data.frame(read_excel("siblings-extended-good.xlsx"))
#################################
profsgood <- profsgood[order(profsgood$refdate),]
seuil <- median(profsgood$titles)
seuil33 <- quantile(profsgood$titles,.33333)
seuil66 <- quantile(profsgood$titles,.66666)
profstop <- profsgood[profsgood$titles>seuil,]
profsbottom <- profsgood[profsgood$titles<=seuil,]
pts<-70 #number of points for which we compute the fertility
sha <- .3  # range of quantiles at each point
#plot of coefficient in recursive estimation
date <- (1:pts)
sign <-(1:pts)
t_test_result <- (1:pts)
p_values <- (1:pts)
t <- (1:pts)
maletop <- (1:pts)
malebottom <- (1:pts)
ns<-nrow(profsgood)
#  splitting at the median
for (i in (1:pts)){
smpl<-profsgood[as.integer(1+(i-1)*ns*(1-sha)/(pts-1)):as.integer((i-1)*ns*(1-sha)/(pts-1)+ns*sha),]
profstop <- smpl[smpl$titles>seuil,]
profsbottom <- smpl[smpl$titles<=seuil,]
maletop[i]=mean(profstop[,"MaleDescendants"],na.rm=T)
malebottom[i]=mean(profsbottom[,"MaleDescendants"],na.rm=T)
date[i]=round(mean(smpl$refdate),0)
t_test_result <- t.test(profstop[,"MaleDescendants"], profsbottom[,"MaleDescendants"])
p_values[i] <- t_test_result$p.value
}
profsgood[!is.na(profsgood$MaleDescendants) & !is.na(profsgood$NumberDescendants) & profsgood$MaleDescendants>profsgood$NumberDescendants,"PKey"]
pdf(file="background_children.pdf", width = 11, height = 7)
plot(t,maletop,lwd=3,type="l",xlab="Year of birth", ylab="fertility (male)",frame.plot=TRUE,ylim=c(1,2.5),
main="Male Children of scholars by level of publications",
sub=paste("Final sample. Nobs=",nrow(profsgood),". ",pts," points, sample=",sha,sep=""),col="coral2",xaxt='n')
axis(1,at=c(1,10,20,30,40,50,60,pts),labels=c(as.integer(date[1]),date[10],date[20],date[30],date[40],date[50],date[60],date[pts]))
lines(t,malebottom,lty=1,col="darkgreen",lwd=3)
legend("topright", legend=c("Above median number of publications", "Below median number of publications"),
col=c("coral2", "darkgreen"), lty=1,lwd=3,
box.lty=0)
#polygon(x=c(t,rev(t)),y=c(high05,rev(low05)), col =  adjustcolor("gray", alpha.f = 0.20), border = NA)
for (i in 1:length(sign)) {
if (p_values[i] > 0.1) {
rect(xleft = t[i] - 0.5, xright = t[i] + 0.5, ybottom = .9, ytop = 2.6, col = adjustcolor("lightpink", alpha.f = 0.3), border = NA, add = TRUE)
} else {
rect(xleft = t[i] - 0.5, xright = t[i] + 0.5, ybottom = .9, ytop = 2.6, col = adjustcolor("lightgreen", alpha.f = 0.3), border = NA, add = TRUE)
}
}
dev.off()
pdf(file="background_children.pdf", width = 11, height = 7)
plot(t,maletop,lwd=3,type="l",xlab="Year of birth", ylab="fertility (male)",frame.plot=TRUE,ylim=c(1,2.5),
main="Male Children of scholars by level of publications",
sub=paste("Final sample. Nobs=",nrow(profsgood),". ",pts," points, sample=",sha,sep=""),col="coral2",xaxt='n')
axis(1,at=c(1,10,20,30,40,50,60,pts),labels=c(as.integer(date[1]),date[10],date[20],date[30],date[40],date[50],date[60],date[pts]))
lines(t,malebottom,lty=4,col="darkgreen",lwd=3)
legend("topright", legend=c("Above median number of publications", "Below median number of publications"),
col=c("coral2", "darkgreen"), lty=c(1,3),lwd=3,
box.lty=0)
#polygon(x=c(t,rev(t)),y=c(high05,rev(low05)), col =  adjustcolor("gray", alpha.f = 0.20), border = NA)
for (i in 1:length(sign)) {
if (p_values[i] <= 0.1) {
rect(xleft = t[i] - 0.5, xright = t[i] + 0.5, ybottom = .9, ytop = 2.6, col = adjustcolor("lightgray", alpha.f = 0.3), border = NA, add = TRUE)
}
}
dev.off()
## remove (almost) everything in the working environment.
rm(list = ls())
# Load RColorBrewer
library(RColorBrewer)
# to import xl files
library(readxl)
#to write excel files
require(writexl)
# Set working directory based on the this.path package, which reliably detects the path of the running script. and projects.
if (!require(this.path)) {
install.packages("this.path")
require(this.path)
}
setwd(dirname(this.path()))
#load data
profs <- as.data.frame(read_excel("siblings-extended.xlsx"))
sink("sink from demographic transition.txt")
#select quality observations
#we drop those with less than 2 sibships and year of death of father is unkwown
profsgood <- profs[profs$SibshipSize>2 | (profs$SibshipSize<3 & !is.na(profs$YearDeathFather)),]
print(paste("criterion 1 leads to drop",nrow(profs)-nrow(profsgood),"observations"))
tp<-nrow(profsgood)
#we drop those with one sibship and one descendent
profsgood <- profsgood[profsgood$SibshipSize>1 | profsgood$NumberDescendants>1,]
print(paste("criterion 2 leads to drop",tp-nrow(profsgood),"observations"))
#table of sibshipsizemale to calibrate phi
table(profsgood$SibshipSizeMale)
summary(profsgood$SibshipSizeMale)
# include a rank variable
profsgood$elder <- ifelse(!is.na(profsgood$RankMale) & profsgood$RankMale==1,1,0)
# include a social class variable
strtop <- gsub("[\r\n]", "","archbishop|knighted|marquess|peer|prince|regent|aristocrat|theologian|linguist|king|title|seigneur|professor|councillor|
bishop|mayor|doctor|rector|general|governor|lord|chief justice|colonel|captain|advocate|baron|physician|baronet|owner|lawyer|author|dr|prosecutor|
lecturer|deputy|magistrate|curator|senator|advisor|king's|medical|scholar|writer|count|advocat|attorney|magister|procurator|commissar|businessman|
physicist|marquis|engineer|manager|politician|landlord|artist|orientalist|entrepreneur|industrialist|duke|geographer|conferenceraad|biskop|anatomist|
consultant|conventcommissary|noble|lieutenant-colonel|confessor|reichsmarschall|comminister|barrister|commandant|councilor|co-director|commonwealth marsh|
prince-episcopal|humanist|magistrat|drtheology|philologist|pedagogue|notaire|banker|exchequer|hetman|geologist|adjutant-general|baroness|lagman|adlad|
eversti|oberst|prepositus|logician|chatelain|cardinal|nobleman|prelate|cartographer|councillor-secretary|statesman|counsellor|vice-admiral|assessor|
chancellor|attorney-postulant|investor|apothecary|judge|architect|antiquary|consul|sir|earl|knight|colonel|surgeon|laird|principal|historian|mathematician|
theologist|admiral|commander|jurist|scholast|notary|dean|marshall|chief marshal|provost|president|chamberlain")
strmid <- gsub("[\r\n]", "","elder|thesaurier|printer|clergyman|inventeur|poet|trumpeter|painter|violinist|oboist|shipmaster|musician|alderman|administrator|
monastic|auditor|reformer|preacher|priest|merchant|pastor|officer|trader|master|superintendent|vicar|secretary|clerk|reverend|minister|lieutenant|teacher|
bailiff|baillif|chaplain|pharmacist|inspector|quartermaster|deacon|manufacturer|skipper|schoolmaster|postmaster|foreman|tutor|councilman|chef|shipowner|bookseller|
treasurer|bourgeois|organist|civil servant|minister|lector|glazier|shipbuilder|surveyor|weighmaster|explorer|syndicus|commissioner|archivist|protocollist|churchman|
controller|turner|librarian|registrar|scribe|greffier|jeweler|wholesaler|burgher|adjutant|marksman|nurse|baccalaureus|papermaker|salesman|cavalry|real estate agent|min. of old kirk|
makler|clothier|kornett|adventurer|orfever|reverend|collector|brigadier|armourer|accounter|churchwarden|sheriff|farmer|hundertman|leaseholder|cornet|counsel")
strbot <- gsub("[\r\n]", "","minstrel|artisan|bierbrouwer|saddler|shoemaker|craftsman|policeman|soldier|cutler|factor|blacksmith|gunsmith|goldsmith|tobacconist|brewer|
cossack|piece caster|shear|cartridge|tailor|fisherman|innkeeper|gardener|baker|grocer|tanner|carpenter|assistant|engraver|maker|cutter|receiver|horseman|
employee|guard|tinsmith|comrade|weaver|coalminer|fishing|miller|slave|schmidt|coppersmith|janitor|beater|churchman|hunting|miner|shopkeeper|caretaker|
shipwright|vigneron|guards|non-commissioned|driver|peasantry|fabrics|burner|military|sweeper|forester|grocery|farm|tractor|lorch|crimond")
strunc <- gsub("[\r\n]", "","syndict|appellate|citizen|ministry|college|chancery|head|collegium|pharmacies|ordynek|consistory|seeker|writter|musices|council|
government|school|trade|administrative|confidant|traveler|buyer|communicator|registar|commercial|advice|secretariat|priech|contractual|senevier|
answerer|receiver|piece caster|shear|cartridge|grocery|farm|aids|crimond|hunting|non-commissioned|
tractor|fabrics|military|lorch")
profsgood$OccupationFather <- tolower(profsgood$OccupationFather)
#occupation dummies
nprof <- nrow(profsgood)
profsgood$soctop <- 0
profsgood$socmid <- 0
profsgood$socbot <- 0
profsgood$socna <- 0
for (i in (1:nprof))
{
if (is.na(profsgood[i,"OccupationFather"])) {profsgood[i,"socna"]<-1
} else {
if (grepl(strtop,profsgood[i,"OccupationFather"])) {profsgood[i,"soctop"]<-1
} else {
if (grepl(strmid,profsgood[i,"OccupationFather"])) {profsgood[i,"socmid"]<-1
} else {
if (grepl(strbot,profsgood[i,"OccupationFather"])) {profsgood[i,"socbot"]<-1
} else {
if (grepl(strunc,profsgood[i,"OccupationFather"])) {profsgood[i,"socna"]<-1}
}}}}
if (profsgood[i,"socna"]+profsgood[i,"soctop"]+profsgood[i,"socmid"]+profsgood[i,"socbot"]==0)
{print(profsgood[i,"OccupationFather"]);
profsgood[i,"socna"]<-1}
}
#check for missing occupation
print("check for missing occupation")
profsgood[profsgood$socna==1 & !is.na(profsgood$OccupationFather),c("PKey","OccupationFather")]
#save the results
write_xlsx(profsgood,"siblings-extended-good.xlsx")
profsgood2 <- profsgood[!is.na(profsgood$MaleDescendants) & !is.na(profsgood$NumberDescendants)
& profsgood$NumberDescendants>0,]
#LET us now do demographics
#parties from wrigley
countwrigley <- c(0.191693856,0.236878826,0.232893927,0.174008955,0.086280695,0.031133555,0.007892533)
# Simple Bar Plot
counts1b <- table(profs$SibshipSizeMale)/nrow(profs)
counts2b <- table(profsgood$SibshipSizeMale)/nrow(profsgood)
pdf("histogram-parities-all.pdf")
barplot(rbind(counts1b[1:7],t(countwrigley)),ylab="Frequency (%)",
xlab="Number of male children",main="Distribution of Parities - males - full sample",beside=T,
legend=c("Profs","English pop (Wrigley)"))
dev.off()
pdf("histogram-parities-good.pdf")
barplot(rbind(counts2b[1:7],t(countwrigley)), ylab="Frequency (%)",
xlab="Number of male children",main="Distribution of Parities - males - final sample",beside=T,
legend=c("Profs","English pop (Wrigley)"))
dev.off()
pdf("frequency-publications.pdf")
hist(profsgood$titles,breaks=c(0,.1,1,2,3,4,5,6,7),freq=T,main="",
ylab="Nb. Profs",xlab="Publications")
abline(v=median(profsgood$titles), lwd=3, lty='dashed')
dev.off()
# Create a migration dummy comparing country of birth and death
profsgood$migrant <- toupper(profsgood$Locations_1_Country2016) == toupper(profsgood$Locations_2_Country2016)
profsgood$sexratio <- profsgood$AmongWhichGirls/profsgood$SibshipSizeMale
# Build a color palette,
coul <- colorRampPalette(brewer.pal(31, "Set3"))(21)
# Pie Chart with Percentages
slices <- c(mean(profs$copenhagen),mean(profs$groningen),mean(profs$alund),mean(profs$ulund),mean(profs$uppsala),
mean(profs$acaddnk),mean(profs$rostock),
mean(profs$kiel),mean(profs$abo),mean(profs$dorpat),mean(profs$amsterdam),
mean(profs$Adanzig),mean(profs$edinburgh),mean(profs$konigsberg),
mean(profs$stockholm),mean(profs$glasgow),mean(profs$greifswald),
mean(profs$Gdanzig),mean(profs$petersburg),mean(profs$franeker),mean(profs$haarlem),
mean(profs$aedinburgh),mean(profs$auppsala),mean(profs$aberdeen),mean(profs$andrews),mean(profs$leiden),
mean(profs$hamburg),mean(profs$trondheim),mean(profs$utrecht),mean(profs$soro))
lbls <- c("U Copenhagen","U Groningen","A Lund","U Lund","U Uppsala","A Danish","U Rostock","U Kiel","U Abo","U Dorpat","C Amsterdam",
"A Danzig","U Edinburgh","U Konigsberg","A Stockholm","U Glasgow","U Greifswald","C Danzig",
"A St-Petersburg","U Franeker","A Haarlem","A Edinburgh","A Uppsala","U Aberdeen","U St Andrews","U Leiden","A Hesse",
"A Trondheim","U Utrecht","G Soro")
pct <- round(slices/sum(slices)*100)
lbls <- paste(lbls, pct) # add percents to labels
lbls <- paste(lbls,"%",sep="") # ad % to labels
pdf(file="share-institutions.pdf", width = 10, height = 7)
pie(slices,labels = lbls, col=coul,main=paste("Pie Chart of Institutions. Nobs=",nrow(profs)))
dev.off()
#pie with urban rural
slices <- c(nrow(profsgood[!is.na(profsgood$Urban) & profsgood$Urban==1,]),
nrow(profsgood[!is.na(profsgood$Urban) & profsgood$Urban==0,]),
nrow(profsgood[is.na(profsgood$Urban),]))
lbls <- c("Urban","Rural","NA")
pct <- round(slices/sum(slices)*100)
lbls <- paste(lbls, pct) # add percents to labels
lbls <- paste(lbls,"%",sep="") # ad % to labels
pdf(file="share-urban-rural.pdf", width = 10, height = 7)
pie(slices,labels = lbls, col=coul,main=paste("Pie Chart of Parental Background. Nobs=",nrow(profs)))
dev.off()
#time series plots
pts<-70 #number of points for which we compute the fertility
sha <- 1-pts/100  # range of quantiles at each point
profsgood2$FemaleDescendants <- profsgood2$NumberDescendants-profsgood2$MaleDescendants
#order by date
profsgood <- profsgood[order(profsgood$refdate),]
profs <- profs[order(profs$refdate),]
results <- data.frame(titlesmed=numeric(),titlesmean=numeric(),longevity=numeric(),fertility=numeric(),fertilitymean=numeric(),
fertilitymean2=numeric(),upperclass=numeric(),middleclass=numeric(),
lowerclass=numeric(),naclass=numeric(),shareurban=numeric(),
sexratio=numeric(),sexratio2=numeric(),sharereject=numeric(),stringsAsFactors = FALSE)
ns<-nrow(profsgood)
nsall<-nrow(profs)
#generate time series
for (i in (1:pts)){
t1 <- profsgood[as.integer(1+(i-1)*ns*(1-sha)/(pts-1)),"refdate"]
t2 <- profsgood[as.integer((i-1)*ns*(1-sha)/(pts-1)+ns*sha),"refdate"]
results[i,"sharereject"] <- 1-
nrow(profsgood[profsgood$refdate>=t1 & profsgood$refdate<t2,])/
nrow(profs[profs$refdate>=t1 & profs$refdate<t2,])
results[i,"shareurban"] <- mean(
profsgood[as.integer(1+(i-1)*ns*(1-sha)/(pts-1)):as.integer((i-1)*ns*(1-sha)/(pts-1)+ns*sha),"Urban"],na.rm=T)
results[i,"titlesmed"] <- median(
profsgood[as.integer(1+(i-1)*ns*(1-sha)/(pts-1)):as.integer((i-1)*ns*(1-sha)/(pts-1)+ns*sha),"titles"])
results[i,"titlesmean"] <- mean(
profsgood[as.integer(1+(i-1)*ns*(1-sha)/(pts-1)):as.integer((i-1)*ns*(1-sha)/(pts-1)+ns*sha),"titles"])
results[i,"upperclass"] <- mean(
profsgood[as.integer(1+(i-1)*ns*(1-sha)/(pts-1)):as.integer((i-1)*ns*(1-sha)/(pts-1)+ns*sha),"soctop"],na.rm=T)
results[i,"middleclass"] <- mean(
profsgood[as.integer(1+(i-1)*ns*(1-sha)/(pts-1)):as.integer((i-1)*ns*(1-sha)/(pts-1)+ns*sha),"socmid"],na.rm=T)
results[i,"lowerclass"] <- mean(
profsgood[as.integer(1+(i-1)*ns*(1-sha)/(pts-1)):as.integer((i-1)*ns*(1-sha)/(pts-1)+ns*sha),"socbot"],na.rm=T)
results[i,"longevity"] <- median(
profsgood[as.integer(1+(i-1)*ns*(1-sha)/(pts-1)):as.integer((i-1)*ns*(1-sha)/(pts-1)+ns*sha),"longevity"],na.rm=T)
results[i,"fertility"] <- 2*median(
profsgood[as.integer(1+(i-1)*ns*(1-sha)/(pts-1)):as.integer((i-1)*ns*(1-sha)/(pts-1)+ns*sha),"SibshipSizeMale"])
results[i,"fertilitymean"] <- 2*mean(
profsgood[as.integer(1+(i-1)*ns*(1-sha)/(pts-1)):as.integer((i-1)*ns*(1-sha)/(pts-1)+ns*sha),"SibshipSizeMale"])
results[i,"fertilitymean2"] <- 2*mean(
profsgood2[as.integer(1+(i-1)*ns*(1-sha)/(pts-1)):as.integer((i-1)*ns*(1-sha)/(pts-1)+ns*sha),"MaleDescendants"])
results[i,"sexratio"] <- sum(profsgood[as.integer(1+(i-1)*ns*(1-sha)/(pts-1)):as.integer((i-1)*ns*(1-sha)/(pts-1)+ns*sha),"AmongWhichGirls"])/
sum(profsgood[as.integer(1+(i-1)*ns*(1-sha)/(pts-1)):as.integer((i-1)*ns*(1-sha)/(pts-1)+ns*sha),"SibshipSizeMale"]-1)
results[i,"sexratio2"] <- sum(profsgood[as.integer(1+(i-1)*ns*(1-sha)/(pts-1)):as.integer((i-1)*ns*(1-sha)/(pts-1)+ns*sha),"FemaleDescendants"])/
sum(profsgood[as.integer(1+(i-1)*ns*(1-sha)/(pts-1)):as.integer((i-1)*ns*(1-sha)/(pts-1)+ns*sha),"MaleDescendants"])
}
results$naclass=1-results$upperclass-results$middleclass-results$lowerclass
date <- (1:pts)
for (i in (1:pts)){
date[i]=round(quantile(profsgood$refdate,i/100+.15)[[1]],0)
}
#export longevity for later treatment
outt <- data.frame(
date = date,
longevity = unlist(results$longevity)
)
write.csv(outt, "UTHC.csv", row.names = FALSE)
pdf(file="socialclass.pdf", width = 10, height = 7)
plot(results[,"upperclass"],type="l",xlab="Year of Birth", ylab="Share",lwd=2,
frame.plot=TRUE,ylim=c(0,.6),col="#000000", lty=1,
main="Proportions of social classes over time",sub=paste("Final sample. Nobs=",nrow(profsgood)),axes=FALSE)
axis(1,at=c(1,10,20,30,40,50,60,pts),labels=c(as.integer(date[1]),date[10],date[20],date[30],date[40],date[50],date[60],date[pts]))
axis(2,at=seq(0,.6, by=.1),labels=c(0,.1,.2,.3,.4,.5,.6))
lines(results[,"middleclass"],col="#0072B2",lwd=2,lty=2)
lines(results[,"lowerclass"],col="#D55E00",lwd=3,lty=3)
lines(results[,"naclass"],col="#009E73",lwd=2,lty=4)
legend(10, .3, legend=c("Upper class", "Middle class","Lower class","NA"),
col=c("#000000","#0072B2","#D55E00","#009E73"), lwd=2, lty=c(1,2,3,4), box.lty=0)
dev.off()
pdf(file="ntitles.pdf", width = 10, height = 7)
plot(results[,"titlesmed"],type="l",xlab="Year of Birth", ylab="Share",lwd=2,col=coul[2],
frame.plot=TRUE,main="Number of titles (asinh)",ylim=c(2.6,4),
sub=paste("Final sample. Nobs=",nrow(profsgood)),axes=FALSE)
axis(1,at=c(1,10,20,30,40,50,60,pts),labels=c(as.integer(date[1]),date[10],date[20],date[30],date[40],date[50],date[60],date[pts]))
axis(2,at=seq(2.6,4, by=.2),labels=c(2.6,2.8,3,3.2,3.4,3.6,3.8,4))
lines(results[,"titlesmean"],col=coul[4],lwd=2)
legend(10,4 , legend=c("Median", "Mean"),
col=c(coul[2],coul[4]),  lwd=2,box.lty=0)
dev.off()
pdf(file="shareurban.pdf", width = 10, height = 7)
plot(results[,"shareurban"],type="l",xlab="Year of Birth", ylab="Share",
frame.plot=TRUE,main="Share of scholars born in cities",ylim=c(.3,.6),
sub=paste("Final sample. Nobs=",nrow(profsgood)),axes=FALSE)
axis(1,at=c(1,10,20,30,40,50,60,pts),labels=c(as.integer(date[1]),date[10],date[20],date[30],date[40],date[50],date[60],date[pts]))
axis(2,at=seq(.3,.6, by=.1),labels=c(0.3,0.4,0.5,0.6))
dev.off()
pdf(file="fertility.pdf", width = 10, height = 7)
plot(results[,"fertility"],type="l",xlab="Year of Birth", ylab="Median sibshipsize",ylim=c(4,8),
frame.plot=TRUE,main="Median size of families",
sub=paste("Final sample. Nobs=",nrow(profsgood)),axes=FALSE)
axis(1,at=c(1,10,20,30,40,50,60,pts),labels=c(as.integer(date[1]),date[10],date[20],date[30],date[40],date[50],date[60],date[pts]))
axis(2,at=seq(5,7, by=1),labels=c(5,6,7))
dev.off()
pdf(file="fertilitymean.pdf", width = 10, height = 7)
plot(results[,"fertilitymean"],type="l",xlab="Year of Birth", ylab="Mean sibshipsize",ylim=c(4,8),
frame.plot=TRUE,main="Mean sibshipsize",
sub=paste("Final sample. Nobs=",nrow(profsgood)),axes=FALSE)
axis(1,at=c(1,10,20,30,40,50,60,pts),labels=c(as.integer(date[1]),date[10],date[20],date[30],date[40],date[50],date[60],date[pts]))
axis(2,at=seq(5,7, by=1),labels=c(5,6,7))
dev.off()
pdf(file="sexratio.pdf", width = 10, height = 7)
plot(results[,"sexratio"],type="l",xlab="Year of Birth", ylab="Sex ratio",
frame.plot=TRUE,ylim=c(.7,1),main="Sex Ratio (F/M)",
sub=paste("Final sample. Nobs=",nrow(profsgood)),axes=FALSE)
axis(1,at=c(1,10,20,30,40,50,60,pts),labels=c(as.integer(date[1]),date[10],date[20],
date[30],date[40],date[50],date[60],date[pts]))
axis(2,at=seq(0.6,1, by=.1),labels=c(.6,.7,.8,.9,1))
dev.off()
tsex_ratio <- rep(0,5)
esex_ratio <- rep(0,5)
proba_boy <- 1.05/2.05
for (i in (1:5))
{tsex_ratio[i] <- (1+proba_boy*(i))/((1-proba_boy)*(i))
esex_ratio[i] <- sum(profsgood[profsgood$SibshipSize==i+1,"SibshipSizeMale"])/sum(profsgood[profsgood$SibshipSize==i+1,"AmongWhichGirls"])
}
print("sex ratio analysis")
print(esex_ratio)
print(tsex_ratio)
#rank analysis
print("rank analysis")
table(profsgood$RankMale)
nrow(profsgood[!is.na(profsgood$RankMale) & profsgood$RankMale==1,])
nrow(profsgood[!is.na(profsgood$RankMale),])
nrow(profsgood[!is.na(profsgood$RankMale) & profsgood$RankMale==1,])/
nrow(profsgood[!is.na(profsgood$RankMale),])
sink()
