######################################################
#
# Bite Force
#
######################################################

## 1. How does the relationship between body mass and bite force in 
## dominant and subordinate naked mole rats compare to the relationship
## seen across the other rodents in the data set?

#### ROSALIND'S ATTEMPT ####

bite <- read.table(file.choose(), header = T, sep = ",") # load bite_force.csv

nmr_dom <- subset(bite, species=="naked mole-rat (dominant)") # subset dominant NMR
nmr_sub <- subset(bite, species=="naked mole-rat (subordinate)") # subset subordinate NMR
rodent <- subset(bite, group == 'rodentia' 
                 & scientific_name !="Heterocephalus glaber") # subset other rodents excluding NMR

summary(nmr_dom$body_mass) # check min and max for axes limits
summary(nmr_dom$bite_force)
summary(nmr_sub$body_mass)
summary(nmr_sub$bite_force)
summary(rodent$body_mass)
summary(rodent$bite_force)

plot(rodent$body_mass,rodent$bite_force, # plot other rodents data
     xlim=c(0,700), ylim=c(0,80), xaxs = "i",yaxs = "i", 
     xlab="Body Mass (g)", ylab="Bite Force (N)", 
     pch=15, col="orange")
lm(rodent$bite_force~rodent$body_mass)
summary(rodent$body_mass)
x0<-6.5
x1<-639
a<-15.67805
b<-0.06377
segments(x0, a+b*x0, x1, a+b*x1, col = "orange", lty = 1, lwd=2) # add linear model line
text(420,10, labels="bite force = 15.68 + (0.064 x mass) ",
     pos=1, cex= 1,font=2,col="orange")
par(new=TRUE) # prepare to add to existing plot
plot(nmr_dom$body_mass,nmr_dom$bite_force, # plot dominant NMR data
     xlim=c(0,700), ylim=c(0,80), xaxs = "i",yaxs = "i", 
     xlab="", ylab="", xaxt="n", yaxt="n",
     pch=21, bg="orange",col="black",cex=1.5,lwd=2)
par(new=TRUE) # prepare to add to existing plot
plot(nmr_sub$body_mass,nmr_sub$bite_force, # plot subordinate NMR data
     xlim=c(0,700), ylim=c(0,80), xaxs = "i",yaxs = "i", 
     xlab="", ylab="", xaxt="n", yaxt="n",
     pch=21, bg="orange",col="darkmagenta",cex=1.5,lwd=2)

text(nmr_dom$bite_force~nmr_dom$body_mass, labels="dom",
     pos=1, cex= 1,font=2,col="black") # add text label to dominant NMR data point
text(nmr_sub$bite_force~nmr_sub$body_mass, labels="sub",
     pos=3, cex= 1,font=2,col="darkmagenta") # add text label to subordinate NMR data point


#### GRAEME'S ATTEMPT ####

bite <- read.table(file.choose(), header = T, sep = ",")

nmr_dom <- subset(bite, species=="naked mole-rat (dominant)")
nmr_sub <- subset(bite, species=="naked mole-rat (subordinate)")
rodent <- subset(bite, group == 'rodentia' & scientific_name !="Heterocephalus glaber")

plot(rodent$body_mass,rodent$bite_force, 
     log = 'x', # plot with a logarithmic x-axis
     xlim=c(5,700), ylim=c(0,80), xaxs = "i",yaxs = "i", 
     xlab="Body Mass (g)", ylab="Bite Force (N)", 
     pch=15, col="green")

xvalues <- seq(from = 5, to = 700, by = 1) # create equally-spaced x-values
yvalues <- 15.68+(xvalues*0.0064) # calculate y-values using linear model equation

par(new = TRUE)
lines(xvalues, yvalues,lty=1,lwd=1,col="green") # plot values as line

par(new=TRUE)
plot(nmr_dom$body_mass,nmr_dom$bite_force, log = 'x',
     xlim=c(5,700), ylim=c(0,80), xaxs = "i",yaxs = "i", 
     xlab="", ylab="", xaxt="n", yaxt="n",
     pch=21, bg="black",col="black",cex=1.5,lwd=2)
par(new=TRUE)
plot(nmr_sub$body_mass,nmr_sub$bite_force, log = 'x',
     xlim=c(5,700), ylim=c(0,80), xaxs = "i",yaxs = "i", 
     xlab="", ylab="", xaxt="n", yaxt="n",
     pch=21, bg="darkmagenta",col="darkmagenta",cex=1.5,lwd=2)

text(nmr_dom$bite_force~nmr_dom$body_mass, labels="dominant",pos=1, cex= 1,font=2,col="black")
text(nmr_sub$bite_force~nmr_sub$body_mass, labels="subordinate",pos=3, cex= 1,font=2,col="darkmagenta")

######################################

## 2. How does the relationship between body mass and bite force in 
## dominant and subordinate naked mole rats compare to the relationship
## seen across the other rodents in the data set?

#### ROSALIND'S ATTEMPT ####

summary(bite$bite_force) # look at bite force variable
summary(bite$body_mass) # look at body mass variable

bite$bite_log<-log(bite$bite_force) # create logged bite force variable
bite$mass_log<-log(bite$body_mass) # create logged body mass variable

nmr_dom <- subset(bite, species=="naked mole-rat (dominant)")
nmr_sub <- subset(bite, species=="naked mole-rat (subordinate)")
rodent <- subset(bite, group == 'rodentia' & scientific_name !="Heterocephalus glaber")
opossum <- subset(bite, group == 'didelphimorphia')
carnivore <- subset(bite, group == 'carnivora')
bat <- subset(bite, group == 'chiroptera')
primate <- subset(bite, group == 'primate') # subset the different groups, 
                                            # including dom and sub NMR

plot(rodent$mass_log,rodent$bite_log, # plot logged data for other rodents
     xlim=c(0,14), ylim=c(0,8), xaxs = "i",yaxs = "i", 
     xlab="Log Body Mass (g)", ylab="Log Max Bite Force (N)", 
     pch=16, col="orange")
lm(rodent$bite_log~rodent$mass_log)
summary(rodent$mass_log)
x0<-1.872
x1<-6.460
a<-0.7703
b<-0.5059
#segments(x0, a+b*x0, x1, a+b*x1, col = "orange", lty = 1, lwd=2)
par(new=TRUE)
plot(opossum$mass_log,opossum$bite_log, # plot logged data for opossums
     xlim=c(0,14), ylim=c(0,8), xaxs = "i",yaxs = "i", 
     xlab="", ylab="", xaxt="n", yaxt="n",
     pch=17, col="hotpink")
lm(opossum$bite_log~opossum$mass_log)
summary(opossum$mass_log)
x0<-4.5
x1<-8.517
a<--0.3681
b<-0.7584
#segments(x0, a+b*x0, x1, a+b*x1, col = "hotpink", lty = 1, lwd=2) 
par(new=TRUE)
plot(carnivore$mass_log,carnivore$bite_log, # plot logged data for carnivores
     xlim=c(0,14), ylim=c(0,8), xaxs = "i",yaxs = "i", 
     xlab="", ylab="", xaxt="n", yaxt="n",
     pch=15, col="royalblue")
lm(carnivore$bite_log~carnivore$mass_log)
summary(carnivore$mass_log)
x0<-6.768
x1<-12.593
a<-0.2644
b<-0.5612
#segments(x0, a+b*x0, x1, a+b*x1, col = "royalblue", lty = 1, lwd=2)  
par(new=TRUE)
plot(bat$mass_log,bat$bite_log, # plot logged data for bats
     xlim=c(0,14), ylim=c(0,8), xaxs = "i",yaxs = "i", 
     xlab="", ylab="", xaxt="n", yaxt="n",
     pch=18, col="forestgreen")
lm(bat$bite_log~bat$mass_log)
summary(bat$mass_log)
x0<-1.386
x1<-7.061
a<--0.9108
b<-0.9030
#segments(x0, a+b*x0, x1, a+b*x1, col = "forestgreen", lty = 1, lwd=2) 
par(new=TRUE)
plot(primate$mass_log,primate$bite_log, # plot logged data for primates
     xlim=c(0,14), ylim=c(0,8), xaxs = "i",yaxs = "i", 
     xlab="", ylab="", xaxt="n", yaxt="n",
     pch=19, col="red")
lm(primate$bite_log~primate$mass_log)
summary(primate$mass_log)
x0<-7.963
x1<-8.081
a<--49.311
b<-7.023
#segments(x0, a+b*x0, x1, a+b*x1, col = "red", lty = 1, lwd=2) 
lm(bite$bite_log~bite$mass_log)
summary(bite$mass_log)
x0<-1.386
x1<-12.593
a<-0.2868
b<-0.5787
segments(x0, a+b*x0, x1, a+b*x1, col = "black", lty = 1, lwd=2)

segments(nmr_dom$mass_log, nmr_dom$bite_log, 
         nmr_dom$mass_log+1, nmr_dom$bite_log+3.25, 
         col = "black", lty = 2, lwd=1) # place label line for NMR dom
text(nmr_dom$mass_log+1, nmr_dom$bite_log+3.75, 
     labels="Dominant NMR",col="black") # place label text for NMR dom
segments(nmr_sub$mass_log, nmr_sub$bite_log, 
         nmr_sub$mass_log-1.5, nmr_sub$bite_log+2.5, 
         col = "darkmagenta", lty = 2, lwd=1) # place label line for NMR sub
text(nmr_sub$mass_log-1.5, nmr_sub$bite_log+3, 
     labels="Subordinate NMR",col="darkmagenta") # place label text for NMR sub

par(new=TRUE)
plot(nmr_dom$mass_log,nmr_dom$bite_log, # plot logged data for NMR dom
     xlim=c(0,14), ylim=c(0,8), xaxs = "i",yaxs = "i", 
     xlab="", ylab="", xaxt="n", yaxt="n",
     pch=21, bg="orange",col="black",cex=1.5,lwd=2)
par(new=TRUE)
plot(nmr_sub$mass_log,nmr_sub$bite_log, # plot logged data for NMR sub
     xlim=c(0,14), ylim=c(0,8), xaxs = "i",yaxs = "i", 
     xlab="", ylab="", xaxt="n", yaxt="n",
     pch=21, bg="orange",col="darkmagenta",cex=1.5,lwd=2)

legend(10,3.5, bg="white", bty="n",pch = c(16,17,15,18,19), 
       col = c("orange", "hotpink","royalblue","forestgreen","red"), 
      legend = c("Rodents","Opossums","Carnivores","Bats","Primates")) # add legend


#### GRAEME'S ATTEMPT ####

bite$bite_log<-log10(bite$bite_force) # create log10 bite force variable
bite$mass_log<-log10(bite$body_mass) # create log10 body mass variable

summary(bite$bite_log)
summary(bite$mass_log)

nmr_dom <- subset(bite, species=="naked mole-rat (dominant)")
nmr_sub <- subset(bite, species=="naked mole-rat (subordinate)")
rodent <- subset(bite, group == 'rodentia' & scientific_name !="Heterocephalus glaber")
opossum <- subset(bite, group == 'didelphimorphia')
carnivore <- subset(bite, group == 'carnivora')
bat <- subset(bite, group == 'chiroptera')
primate <- subset(bite, group == 'primate')

plot(rodent$mass_log,rodent$bite_log, 
     xlim=c(0,4), ylim=c(0,3), xaxs = "i",yaxs = "i", 
     xlab="Log10 Body Mass (g)", ylab="Log10 Max Bite Force (N)", 
     pch=16, col="orange")
par(new=TRUE)
plot(opossum$mass_log,opossum$bite_log, 
     xlim=c(0,4), ylim=c(0,3), xaxs = "i",yaxs = "i", 
     xlab="", ylab="", xaxt="n", yaxt="n",
     pch=17, col="hotpink")
par(new=TRUE)
plot(carnivore$mass_log,carnivore$bite_log, 
     xlim=c(0,4), ylim=c(0,3), xaxs = "i",yaxs = "i", 
     xlab="", ylab="", xaxt="n", yaxt="n",
     pch=15, col="royalblue")
par(new=TRUE)
plot(bat$mass_log,bat$bite_log, 
     xlim=c(0,4), ylim=c(0,3), xaxs = "i",yaxs = "i", 
     xlab="", ylab="", xaxt="n", yaxt="n",
     pch=18, col="forestgreen")
par(new=TRUE)
plot(primate$mass_log,primate$bite_log, 
     xlim=c(0,4), ylim=c(0,3), xaxs = "i",yaxs = "i", 
     xlab="", ylab="", xaxt="n", yaxt="n",
     pch=19, col="red")
lm(bite$bite_log~bite$mass_log)
summary(bite$mass_log)
x0<-0.602
x1<-5.47
a<-0.1246
b<-0.5787
segments(x0, a+b*x0, x1, a+b*x1, col = "black", lty = 1, lwd=2)

segments(nmr_dom$mass_log, nmr_dom$bite_log, 
         nmr_dom$mass_log+0.2, nmr_dom$bite_log+1, 
         col = "black", lty = 2, lwd=1)
text(nmr_dom$mass_log+0.2, nmr_dom$bite_log+1, 
     labels="Dominant NMR",col="black")
segments(nmr_sub$mass_log, nmr_sub$bite_log, 
         nmr_sub$mass_log-0.5, nmr_sub$bite_log+1.5, 
         col = "darkmagenta", lty = 2, lwd=1)
text(nmr_sub$mass_log-0.5, nmr_sub$bite_log+1.5, 
     labels="Subordinate NMR",col="darkmagenta")

par(new=TRUE)
plot(nmr_dom$mass_log,nmr_dom$bite_log, 
     xlim=c(0,4), ylim=c(0,3), xaxs = "i",yaxs = "i", 
     xlab="", ylab="", xaxt="n", yaxt="n",
     pch=21, bg="orange",col="black",cex=1.5,lwd=2)
par(new=TRUE)
plot(nmr_sub$mass_log,nmr_sub$bite_log, 
     xlim=c(0,4), ylim=c(0,3), xaxs = "i",yaxs = "i", 
     xlab="", ylab="", xaxt="n", yaxt="n",
     pch=21, bg="orange",col="darkmagenta",cex=1.5,lwd=2)

legend(2.5,0.9, bg="white", bty="n",pch = c(16,17,15,18,19),
       col = c("orange", "hotpink","royalblue","forestgreen","red"), 
       legend=c("Rodents","Opossums","Carnivores","Bats","Primates"))


######################################################
#
# Women in STEM
#
######################################################

#### ROSALIND'S ATTEMPT ####

stem <- read.table(file.choose(), header = T, sep = ",") # load core_stem.csv
View(stem) # view the data

summary(stem$year) # check min and max of variables
summary(stem$women_percent)
summary(stem$men_percent)

plot(stem$year,stem$women_percent, # plot the women's % over time
     xlim=c(2008,2020), ylim=c(0,100),
     xaxs = "i",yaxs = "i", xlab="Year", ylab="% UK Core STEM workforce", 
     pch=16, col="red")
par(new=TRUE) # prepare to add to existing plot
plot(stem$year,stem$men_percent, # plot the men's % over time
     xlim=c(2008,2020), ylim=c(0,100),
     xaxs = "i",yaxs = "i", xlab="", ylab="", xaxt="n",yaxt="n",
     pch=17, col="seagreen1")
sortedstem <- stem[order(stem$year),] # sort data rows by year
lines(sortedstem$year,sortedstem$women_percent,
      lty=1,lwd=1,col="red") # connect women data points over time
lines(sortedstem$year,sortedstem$men_percent,
      lty=1,lwd=1,col="seagreen1") # connect men data points over time

text(sortedstem$women_percent~sortedstem$year, 
     labels=sortedstem$women_millions, # add total counts of women as labels
     pos=1, cex= 0.8,font=1,col="darkred")
text(sortedstem$men_percent~sortedstem$year, 
     labels=sortedstem$men_millions, # add total counts of men as labels
     pos=3, cex= 0.8,font=1,col="forestgreen")

legend(2009,65, bg="white", bty="n",pch = c(17,16),
       col = c("seagreen1", "red"), 
       legend=c("Men","Women"),lty=1) # add a legend
text(2019.3,88, labels="*",pos=1, cex= 1,col="blue")
text(2019.2,25, labels="*",pos=1, cex= 1,col="blue")
text(2010.1,45, labels="* total count (millions)",
     cex= 0.8,col="blue") # add text to explain the total count labels


#### GRAEME'S ATTEMPT ####

stem <- read.table(file.choose(), header = T, sep = ",")

plot(stem$year,stem$women_percent, # just plotting women data
     xlim=c(2008,2020), 
     ylim=c(15.8,25.2), # reduce y-axis range
     xaxs = "i",yaxs = "i", xlab="Year", ylab="% Women in STEM", 
     pch=16, col="red")

sortedstem <- stem[order(stem$year),]
lines(sortedstem$year,sortedstem$women_percent,lty=1,lwd=1,col="red")

abline(h=(seq(16,25,0.25)), col="lightgray", lty=1) # add minor grid lines
abline(h=(seq(16,25,1)), col="darkgray", lty=1, lwd = 2) # add major grid lines

text(sortedstem$women_percent~sortedstem$year, 
     labels=sortedstem$women_millions, # add total counts of women as labels
     pos=1, cex= 0.8,font=1,col="darkred")