facialhairstyles.R

library(gsheet)
library(tidyverse)
library(magrittr)
library(fANCOVA)
library(Cairo)
#library(plotly)
facialhairstyles <-
  gsheet::gsheet2tbl(
    "https://docs.google.com/spreadsheets/d/14hgOhTHb5caVmmLwq859U22BlsIBeQyGPV1yAJ74Xms/edit#gid=0"
  )
facialhairstyles %<>% transmute(
  Year = X1,
  Cleanshaven,
  Sideburns,
  Muttonchops,
  Moustache,
  `Moustache (with Sideburns)`,
  `Moustache (with Muttonchops)`,
  Goatee,
  `Full Beard/Chinstrap`
)

# Brillouin's index of evenness
Hb2 <- function(ns) {
  N <- sum(ns)
  S <- length(ns)
  n <- floor(N/S)
  r <- N-S*n
  ((lfactorial(N) - sum(lfactorial(ns)))/(lfactorial(N)-r*lfactorial(n+1)-(S-r)*lfactorial(n)))
}

diversity <- 1:length(facialhairstyles$Year)

for (i in 1:length(facialhairstyles$Year)) {
  diversity[i] <- Hb2(facialhairstyles[i, 2:9])
}
diversity <-
  cbind(facialhairstyles$Year, diversity) %>% as.data.frame()
names(diversity)[1] <- "year"

test <-
  loess.as(x = diversity$year,
           y = diversity$diversity,
           degree = 2)
test$pars$span
plotobj <-
  diversity %>%
  ggplot(aes(x = year,
             y = diversity)) +
  geom_line() +
  labs(x = "Year",
       y = "Brillouin's standardized index of diversity (evenness)") +
  ggtitle("Diversity of University Graduate Facial Hair Styles 1898-2008")+
  geom_smooth(span=loess.as(x = diversity$year,
                            y = diversity$diversity,
                            degree = 2)$pars$span)+
  scale_y_continuous(limits=c(-0.01,1))+
  scale_x_continuous(breaks=seq(1900,2010,10))
Cairo()
plotobj
dev.off()
#ggplotly(plotobj)