
## 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()