Zebra Teeth Simulation Outcomes

Author

Ben Davies, Alex Norwood, Tyler Faith

Load packages

library(tidyverse)

Birth and mortality profiles

Birth

birthDeath<-read_csv("birthDeath.csv")
birthDeath$born<-month.abb[birthDeath$born]

ggplot(birthDeath,aes(x=fct_relevel(as_factor(born),"Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"))) +
  geom_bar(fill="darkblue") +
  theme_classic() +
  labs(x="Birth month",y="n")  +
  theme(
    axis.text=element_text(size=12),
    axis.title=element_text(size=14),
    axis.text.x = element_text(angle = 45,hjust=1))

Death

ggplot(birthDeath,aes(x=died)) +
  geom_histogram(fill="darkred") +
  theme_classic() +
  labs(x="Age at death (years)",y="n") +
  theme(
    axis.text=element_text(size=12),
    axis.title=element_text(size=14))

Tooth growth models

toothGrowth<-read_csv("growthRates.csv") %>%
  pivot_longer(cols=constant:nonLinear,names_to="Model",values_to="crownHeight")
ggplot(toothGrowth,aes(x=day,y=crownHeight,color=fct_recode(Model,Constant="constant","Non-linear"="nonLinear"))) +
  geom_line(lwd=1.1) +
  scale_color_manual(values=c("black","darkorange")) +
  theme_classic() +
  labs(x="Days since start of mineralization",y="Crown Height (\u{03BC}m)",color="Model") +
  theme(axis.text=element_text(size=12),axis.title=element_text(size=14))

Model 𝝳¹⁸O signals

Amplitudes of 1, 2, 3, and 4‰ were simulated; most results below are based on the 3‰ signal.

d18O<-read_csv("sim_d18O.csv") %>%
  pivot_longer(cols=AMP1:AMP4,names_to="Amplitude",values_to="d18O")
ggplot(d18O,aes(x=as_date(DAY,origin="2000-12-21"),y=d18O,color=case_match(Amplitude,"AMP1"~"1","AMP2"~"2","AMP3"~"3","AMP4"~"4"))) +
  geom_hline(yintercept=-4,linetype=2)+
  geom_line(lwd=1.1) +
  scale_color_brewer(palette="Set1") +
  scale_x_date(date_labels = "%b", date_breaks = '3 months') +
  theme_classic() +
  labs(x="",y=expression(paste(delta^"18","O")),color="Amplitude (\u{2030})") +
  theme(axis.text=element_text(size=12),axis.title=element_text(size=14))

Constant vs non-linear growth rate (no wear; Amplitude = 3)

These teeth come from animals that are all born on the same day (Jan 1) and die after 55 months.

Just the constant model:

Ex1<-read_csv("simZebraDataEX12.csv") %>% pivot_longer(cols=d18O_1:d18O_4,names_to="Amplitude",values_to="d18O",names_prefix="d18O_")
Ex1_Amp3<-filter(Ex1,Amplitude==3,specimen==1)
Ex1_Amp3_con<-filter(Ex1_Amp3,model=="Constant")

ggplot(Ex1_Amp3_con,aes(x=mm_band,y=d18O)) + 
  geom_line() +
   geom_hline(yintercept=-4,linetype=2)+
  geom_line(lwd=1.1) +
  theme_classic() +
  labs(x="Distance from Occlusal Surface (mm)",y=expression(paste(delta^"18","O")))+
  theme(axis.text=element_text(size=12),axis.title=element_text(size=14))

Constant and non-linear

ggplot(Ex1_Amp3,aes(x=mm_band,y=d18O,color=fct_recode(model,"Non-Linear" = "High"))) + geom_line() +
   geom_hline(yintercept=-4,linetype=2)+
  geom_line(lwd=1.1) +
  scale_color_manual(values=c("black","darkorange")) +
  theme_classic() +
  labs(x="Distance from Occlusal Surface (mm)",y=expression(paste(delta^"18","O")),color="Model")+
  theme(axis.text=element_text(size=12),axis.title=element_text(size=14))

Non-linear growth with variable birth season ( Amplitude = 3)

Ex3<-read_csv("simZebraDataEX3.csv") %>% pivot_longer(cols=d18O_1:d18O_4,names_to="Amplitude",values_to="d18O",names_prefix="d18O_")

Rows: 55000 Columns: 17
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (1): model
dbl (16): specimen, age_at_death, birth_year, birth_month, birth_day, crown_...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

uv<-unique(Ex3$specimen)
samp<-sample(uv,50,replace=FALSE)
Ex3_Amp3<-filter(Ex3,Amplitude==3,model=="High", specimen %in% samp)

ggplot(Ex3_Amp3,aes(x=mm_band,y=d18O,group=specimen,color=specimen)) + 
   geom_hline(yintercept=-4,linetype=2)+
  geom_line(alpha=0.8,color="grey") +
  theme_classic() +
  labs(x="Distance from Occlusal Surface (mm)",y=expression(paste(delta^"18","O"))) +
  theme(axis.text=element_text(size=12),axis.title=element_text(size=14))

Non-linear growth with variable birth season and death age (Amplitude = 3)

Note: samples chosen random

set.seed(303) #to keep the sampling of teeth consistent
Ex4<-read_csv("simZebraDataEX4.csv")

Rows: 33671 Columns: 17
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (1): model
dbl (16): specimen, age_at_death, birth_year, birth_month, birth_day, crown_...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Ex4_pivot<-Ex4 %>% pivot_longer(cols=d18O_1:d18O_4,names_to="Amplitude",values_to="d18O",names_prefix="d18O_")
uv<-unique(Ex4$specimen)
samp<-sample(uv,5,replace=FALSE)
Ex4_Amp3_samp<-filter(Ex4_pivot,Amplitude==3,model=="High", specimen %in% samp)

ggplot(Ex4_Amp3_samp,aes(x=mm_band,y=d18O,group=as_factor(specimen),color=as_factor(specimen))) + 
   geom_hline(yintercept=-4,linetype=2)+
  geom_line(lwd=1.1) +
  theme_classic() +
  scale_color_brewer(palette="Set1") +
  labs(x="Distance from Occlusal Surface (mm)",y=expression(paste(delta^"18","O")),color="Specimen #")+
  theme(axis.text=element_text(size=12),axis.title=element_text(size=14))

Boxplots if intratooth variation

ggplot(Ex4_Amp3_samp,aes(x=as_factor(specimen),y=d18O,group=as_factor(specimen),color=as_factor(specimen))) + 
   geom_hline(yintercept=-4,linetype=2)+
  geom_boxplot(show.legend = FALSE) +
  geom_hline(yintercept=-4,linetype=2)+
  theme_classic() +
  scale_color_brewer(palette="Set1") +
  labs(x="Specimen #",y=expression(paste(delta^"18","O")),color="Specimen #")+
  theme(axis.text=element_text(size=12),axis.title=element_text(size=14))

BULK SAMPLING METHOD A: Estimating amplitude from random 1mm bands (Non-linear growth, Amplitude = 3)

Ex4_high<-Ex4 %>%  filter(model=="High") %>% drop_na()

sampleSizes<-c(10,20,50,100)
size<-c()
bootstrapAmp<-c()
for (i in sampleSizes) {
  for (j in c(1:100)) {
  samp<-Ex4_high %>% group_by(specimen) %>% slice_sample(n=1) %>% ungroup() %>% slice_sample(n=i)
  size<-append(size,i)
  bootstrapAmp<-append(bootstrapAmp,abs(diff(range(samp$d18O_3)))/2)
  }
}
data<-data.frame(size,bootstrapAmp)
ggplot(data,aes(x=as_factor(size),y=bootstrapAmp)) +
  geom_hline(yintercept=3,linetype=2)+
  geom_boxplot() +
  theme_classic() +
  ylim(c(1,4)) +
  labs(x="Sample Size (# of teeth)",y="Estimated Amplitude (\u{2030})") +
  theme(axis.text=element_text(size=12),axis.title=element_text(size=14))

BULK SAMPLING METHOD A: Estimating amplitude from sampling a random 1mm band with 20mm size cutoff (Non-linear growth, Amplitude = 3)

cutoff<-20
Ex4_highCut<-Ex4 %>%  filter(model=="High") %>% drop_na() %>% group_by(specimen) %>% filter( n() >=cutoff) %>% ungroup()

sampleSizes<-c(10,20,50,100)
size<-c()
bootstrapAmp<-c()
for (i in sampleSizes) {
  for (j in c(1:100)) {
  samp<-Ex4_highCut %>% group_by(specimen) %>% slice_sample(n=1) %>% ungroup() %>% slice_sample(n=i)
  size<-append(size,i)
  bootstrapAmp<-append(bootstrapAmp,abs(diff(range(samp$d18O_3)))/2)
  }
}
data<-data.frame(size,bootstrapAmp)
ggplot(data,aes(x=as_factor(size),y=bootstrapAmp)) +
  geom_hline(yintercept=3,linetype=2)+
  geom_boxplot() +
  theme_classic() +
  ylim(c(1,4)) +
  labs(x="Sample Size (# of teeth)",y="Estimated Amplitude (\u{2030})") +
  theme(axis.text=element_text(size=12),axis.title=element_text(size=14))

BULK SAMPLING METHOD B: Estimating amplitude from averaged sample taken along growth axis (Non-linear growth, Amplitude = 3)

Ex4_high<-Ex4 %>%  filter(model=="High") %>% drop_na()

sampleSizes<-c(10,20,50,100)
size<-c()
bootstrapAmp<-c()
for (i in sampleSizes) {
  for (j in c(1:100)) {
  samp<-Ex4_high %>% group_by(specimen) %>% slice_sample(n=i) %>%  summarize(mean(d18O_3))
  size<-append(size,i)
  bootstrapAmp<-append(bootstrapAmp,abs(diff(range(samp$`mean(d18O_3)`)))/2)
  }
}
data<-data.frame(size,bootstrapAmp)
ggplot(data,aes(x=as_factor(size),y=bootstrapAmp)) +
  geom_hline(yintercept=3,linetype=2)+
  geom_boxplot() +
  theme_classic() +
  ylim(c(1,4)) +
  labs(x="Sample Size (# of teeth)",y="Estimated Amplitude (\u{2030})") +
  theme(axis.text=element_text(size=12),axis.title=element_text(size=14))

BULK SAMPLING METHOD B: Estimating amplitude from averaged sample taken along growth axis with 20 mm cutoff (Non-linear growth, Amplitude = 3)

cutoff<-20
Ex4_highCut<-Ex4 %>%  filter(model=="High") %>% drop_na() %>% group_by(specimen) %>% filter( n() >=cutoff) %>% ungroup() 

sampleSizes<-c(10,20,50,100)
size<-c()
bootstrapAmp<-c()
for (i in sampleSizes) {
  for (j in c(1:100)) {
  samp<-Ex4_highCut %>% group_by(specimen) %>% slice_sample(n=i) %>%  summarize(mean(d18O_3))
  size<-append(size,i)
  bootstrapAmp<-append(bootstrapAmp,abs(diff(range(samp$`mean(d18O_3)`)))/2)
  }
}
data<-data.frame(size,bootstrapAmp)
ggplot(data,aes(x=as_factor(size),y=bootstrapAmp)) +
  geom_hline(yintercept=3,linetype=2)+
  geom_boxplot() +
  theme_classic() +
  ylim(c(1,4)) +
  labs(x="Sample Size (# of teeth)",y="Estimated Amplitude (\u{2030})") +
  theme(axis.text=element_text(size=12),axis.title=element_text(size=14))