###Exaeretodon allometry modeling
Exa_data <- read.csv("Exaeretodon measurements_log10_09-06-2021.csv")

library(lme4)
library(lmerTest)
library(afex)
library(pbkrtest)
library(plotrix)

#please note that the appropriate spreadsheet to be used here contains log-transformed data, and not the raw mm measurements in the other file

###Checking if data are normally distributed----
with(Exa_data, shapiro.test(BSL)) 
with(Exa_data, shapiro.test(MUL)) 
with(Exa_data, shapiro.test(PAL)) 
with(Exa_data, shapiro.test(OL)) 
with(Exa_data, shapiro.test(IO)) 
with(Exa_data, shapiro.test(TEL))
with(Exa_data, shapiro.test(SW))
with(Exa_data, shapiro.test(BW))
with(Exa_data, shapiro.test(UP))
with(Exa_data, shapiro.test(PD))
with(Exa_data, shapiro.test(TP))
with(Exa_data, shapiro.test(OD))
with(Exa_data, shapiro.test(OW))
with(Exa_data, shapiro.test(BB))
with(Exa_data, shapiro.test(ZH))
with(Exa_data, shapiro.test(DL))
with(Exa_data, shapiro.test(ZW))

####Muzzle----
MUL.u <- Exa_data[!(Exa_data$MUL.D==1),]
lm_basic_MUL <- lm(MUL ~ BSL, data=MUL.u, offset=1.00*BSL)

I <- rep(0, length(Exa_data$MUL), na.rm=TRUE)
I[Exa_data$MUL.D==1] <- 1:sum((Exa_data$MUL.D), na.rm=TRUE)
lmer_basic_MUL <- lmer(Exa_data$MUL ~ Exa_data$BSL + (1|I), offset=1.00*Exa_data$BSL)

summary(lmer_basic_MUL)
summary(lm_basic_MUL)

x1 <- seq(2.0, 2.7, 0.05)

y1.mul <- (1+summary(lmer_basic_MUL)$coef[2,1])*x1 + summary(lmer_basic_MUL)$coef[1,1]

lmer_basic_MUL <- lmer(Exa_data$MUL ~ Exa_data$BSL + (1|I))
plot(Exa_data$MUL ~ Exa_data$BSL, pch=16, cex=1, col="Red", 
     xlab="Skull Length (Log[mm])", ylab="Muzzle Length (Log[mm])")
newx.mul <- seq(min(MUL.u$BSL), max(MUL.u$BSL), length.out=9)
abline(lm(MUL.u$MUL ~ MUL.u$BSL), col="black", lw=2)
conf_interval <- predict(lm_basic_MUL, newdata.mul=data.frame(x=newx.mul), interval="confidence", level=0.95)
x.mul <- MUL.u$BSL
lines(x.mul[order(x.mul)], conf_interval[order(x.mul),2], lty=2, col="Black", lw=2)
lines(x.mul[order(x.mul)], conf_interval[order(x.mul),3], lty=2, col="Black", lw=2)
lines(y1.mul ~ x1, col="Blue", lw=2)
m.mul <- lmer_basic_MUL
newdat.mul <- data.frame(x=seq(2.2, 2.7, length.out = 23))
mm.mul <- model.matrix(~x, newdat.mul)
newdat.mul$y <- mm.mul%*%fixef(m.mul)
pvar1.mul <- diag(mm.mul %*% tcrossprod(vcov(m.mul),mm.mul))
tvar1.mul <- pvar1.mul+VarCorr(m.mul)$f[1]
newdat.mul <- data.frame(
        newdat.mul
        , plo.mul = newdat.mul$y-1.96*sqrt(pvar1.mul)
        , phi.mul = newdat.mul$y+1.96*sqrt(pvar1.mul)
)
lines(newdat.mul$x, newdat.mul$plo.mul, col="Blue", lty=2, lw=2)
lines(newdat.mul$x, newdat.mul$phi.mul, col="Blue", lty=2, lw=2)
points(MUL.u$MUL ~ MUL.u$BSL, col="black", pch=16)

###Palate length----
PAL.u <- Exa_data[!(Exa_data$PAL.D==1),]
lm_basic_PAL <- lm(PAL ~ BSL, data=PAL.u, offset=1.00*BSL)

I <- rep(0, length(Exa_data$PAL), na.rm=TRUE)
I[Exa_data$PAL.D==1] <- 1:sum((Exa_data$PAL.D), na.rm=TRUE)
lmer_basic_PAL <- lmer(Exa_data$PAL ~ Exa_data$BSL + (1|I), offset=1.00*Exa_data$BSL)

summary(lmer_basic_PAL)
summary(lm_basic_PAL)

y1.io <- (1+summary(lmer_basic_PAL)$coef[2,1])*x1 + summary(lmer_basic_PAL)$coef[1,1]

lmer_basic_PAL <- lmer(Exa_data$PAL ~ Exa_data$BSL + (1|I))
plot(Exa_data$PAL ~ Exa_data$BSL, pch=16, cex=1, col="Red", 
     xlab="Skull Length (Log[mm])", ylab="Palate Length (Log[mm])")
newx.pal <- seq(min(PAL.u$BSL), max(PAL.u$BSL), length.out=4)
abline(lm(PAL.u$PAL ~ PAL.u$BSL), col="black", lw=2)
conf_interval <- predict(lm_basic_PAL, newdata.pal=data.frame(x=newx.pal), interval="confidence", level=0.95)
x.pal <- PAL.u$BSL
lines(x.pal[order(x.pal)], conf_interval[order(x.pal),2], lty=2, col="Black", lw=2)
lines(x.pal[order(x.pal)], conf_interval[order(x.pal),3], lty=2, col="Black", lw=2)
lines(y1.pal ~ x1, col="Blue", lw=2)
m.pal <- lmer_basic_PAL
newdat.pal <- data.frame(x=seq(2.2, 2.7, length.out = 23))
mm.pal <- model.matrix(~x, newdat.pal)
newdat.pal$y <- mm.pal%*%fixef(m.pal)
pvar1.pal <- diag(mm.pal %*% tcrossprod(vcov(m.pal),mm.pal))
tvar1.pal <- pvar1.pal+VarCorr(m.pal)$f[1]
newdat.pal <- data.frame(
        newdat.pal
        , plo.pal = newdat.pal$y-1.96*sqrt(pvar1.pal)
        , phi.pal = newdat.pal$y+1.96*sqrt(pvar1.pal)
)
lines(newdat.pal$x, newdat.pal$plo.pal, col="Blue", lty=2, lw=2)
lines(newdat.pal$x, newdat.pal$phi.pal, col="Blue", lty=2, lw=2)
points(PAL.u$PAL ~ PAL.u$BSL, col="black", pch=16)

###Orbit length----
OL.u <- Exa_data[!(Exa_data$OL.D==1),]
lm_basic_OL <- lm(OL ~ BSL, data=OL.u, offset=1.00*BSL)

I <- rep(0, length(Exa_data$OL))
I[Exa_data$OL.D==1] <- 1:sum(Exa_data$OL.D)
lmer_basic_OL <- lmer(Exa_data$OL ~ Exa_data$BSL + (1|I), offset=1.00*Exa_data$BSL)

summary(lmer_basic_OL)
summary(lm_basic_OL)

y1.ol <- (1+summary(lmer_basic_OL)$coef[2,1])*x1 + summary(lmer_basic_OL)$coef[1,1]

lmer_basic_OL <- lmer(Exa_data$OL ~ Exa_data$BSL + (1|I))
plot(Exa_data$OL ~ Exa_data$BSL, pch=16, cex=1, col="Red", 
     xlab="Skull Length (Log[mm])", ylab="Orbit Length (Log[mm])")
newx.ol <- seq(min(OL.u$BSL), max(OL.u$BSL), length.out=4)
abline(lm(OL.u$OL ~ OL.u$BSL), col="black", lw=2)
conf_interval <- predict(lm_basic_OL, newdata.ol=data.frame(x=newx.ol), interval="confidence", level=0.95)
x.ol <- OL.u$BSL
lines(x.ol[order(x.ol)], conf_interval[order(x.ol),2], lty=2, col="Black", lw=2)
lines(x.ol[order(x.ol)], conf_interval[order(x.ol),3], lty=2, col="Black", lw=2)
lines(y1.ol ~ x1, col="Blue", lw=2)
m.ol <- lmer_basic_OL
newdat.ol <- data.frame(x=seq(2.2, 2.7, length.out = 23))
mm.ol <- model.matrix(~x, newdat.ol)
newdat.ol$y <- mm.ol%*%fixef(m.ol)
pvar1.ol <- diag(mm.ol %*% tcrossprod(vcov(m.ol),mm.ol))
tvar1.ol <- pvar1.ol+VarCorr(m.ol)$f[1]
newdat.ol <- data.frame(
        newdat.ol
        , plo.ol = newdat.ol$y-1.96*sqrt(pvar1.ol)
        , phi.ol = newdat.ol$y+1.96*sqrt(pvar1.ol)
)
lines(newdat.ol$x, newdat.ol$plo.ol, col="Blue", lty=2, lw=2)
lines(newdat.ol$x, newdat.ol$phi.ol, col="Blue", lty=2, lw=2)
points(OL.u$OL ~ OL.u$BSL, col="black", pch=16)

###Interorbital distance----
IO.u <- Exa_data[!(Exa_data$IO.D==1),]
lm_basic_IO <- lm(IO ~ BSL, data=IO.u, offset=1.00*BSL)

I <- rep(0, length(Exa_data$IO))
I[Exa_data$IO.D==1] <- 1:sum(Exa_data$IO.D)
lmer_basic_IO <- lmer(Exa_data$IO ~ Exa_data$BSL + (1|I), offset=1.00*Exa_data$BSL)

summary(lmer_basic_IO)
summary(lm_basic_IO)

y1.io <- (1+summary(lmer_basic_IO)$coef[2,1])*x1 + summary(lmer_basic_IO)$coef[1,1]

lmer_basic_IO <- lmer(Exa_data$IO ~ Exa_data$BSL + (1|I))
plot(Exa_data$IO ~ Exa_data$BSL, pch=16, cex=1, col="Red", 
     xlab="Skull Length (Log[mm])", ylab="Interorbital distance (Log[mm])")
newx.io <- seq(min(IO.u$BSL), max(IO.u$BSL), length.out=4)
abline(lm(IO.u$IO ~ IO.u$BSL), col="black", lw=2)
conf_interval <- predict(lm_basic_IO, newdata.io=data.frame(x=newx.io), interval="confidence", level=0.95)
x.io <- IO.u$BSL
lines(x.io[order(x.io)], conf_interval[order(x.io),2], lty=2, col="Black", lw=2)
lines(x.io[order(x.io)], conf_interval[order(x.io),3], lty=2, col="Black", lw=2)
lines(y1.io ~ x1, col="Blue", lw=2)
m.io <- lmer_basic_IO
newdat.io <- data.frame(x=seq(2.2, 2.7, length.out = 23))
mm.io <- model.matrix(~x, newdat.io)
newdat.io$y <- mm.io%*%fixef(m.io)
pvar1.io <- diag(mm.io %*% tcrossprod(vcov(m.io),mm.io))
tvar1.io <- pvar1.io+VarCorr(m.io)$f[1]
newdat.io <- data.frame(
        newdat.io
        , plo.io = newdat.io$y-1.96*sqrt(pvar1.io)
        , phi.io = newdat.io$y+1.96*sqrt(pvar1.io)
)
lines(newdat.io$x, newdat.io$plo.io, col="Blue", lty=2, lw=2)
lines(newdat.io$x, newdat.io$phi.io, col="Blue", lty=2, lw=2)
points(IO.u$IO ~ IO.u$BSL, col="black", pch=16)

###Temporal length----
TEL.u <- Exa_data[!(Exa_data$TEL.D==1),]
lm_basic_TEL <- lm(TEL ~ BSL, data=TEL.u, offset=1.00*BSL)

I <- rep(0, length(Exa_data$TEL))
I[Exa_data$TEL.D==1] <- 1:sum(Exa_data$TEL.D)
lmer_basic_TEL <- lmer(Exa_data$TEL ~ Exa_data$BSL + (1|I), offset=1.00*Exa_data$BSL)

summary(lmer_basic_TEL)
summary(lm_basic_TEL)

y1.tel <- (1+summary(lmer_basic_TEL)$coef[2,1])*x1 + summary(lmer_basic_TEL)$coef[1,1]

lmer_basic_TEL <- lmer(Exa_data$TEL ~ Exa_data$BSL + (1|I))
plot(Exa_data$TEL ~ Exa_data$BSL, pch=16, cex=1, col="Red", 
     xlab="Skull Length (Log[mm])", ylab="Temporal Length (Log[mm])")
newx.tel <- seq(min(TEL.u$BSL), max(TEL.u$BSL), length.out=4)
abline(lm(TEL.u$TEL ~ TEL.u$BSL), col="black", lw=2)
conf_interval <- predict(lm_basic_TEL, newdata.tel=data.frame(x=newx.tel), interval="confidence", level=0.95)
x.tel <- TEL.u$BSL
lines(x.tel[order(x.tel)], conf_interval[order(x.tel),2], lty=2, col="Black", lw=2)
lines(x.tel[order(x.tel)], conf_interval[order(x.tel),3], lty=2, col="Black", lw=2)
lines(y1.tel ~ x1, col="Blue", lw=2)
m.tel <- lmer_basic_TEL
newdat.tel <- data.frame(x=seq(2.2, 2.7, length.out = 23))
mm.tel <- model.matrix(~x, newdat.tel)
newdat.tel$y <- mm.tel%*%fixef(m.tel)
pvar1.tel <- diag(mm.tel %*% tcrossprod(vcov(m.tel),mm.tel))
tvar1.tel <- pvar1.tel+VarCorr(m.tel)$f[1]
newdat.tel <- data.frame(
        newdat.tel
        , plo.tel = newdat.tel$y-1.96*sqrt(pvar1.tel)
        , phi.tel = newdat.tel$y+1.96*sqrt(pvar1.tel)
)
lines(newdat.tel$x, newdat.tel$plo.tel, col="Blue", lty=2, lw=2)
lines(newdat.tel$x, newdat.tel$phi.tel, col="Blue", lty=2, lw=2)
points(TEL.u$TEL ~ TEL.u$BSL, col="black", pch=16)

###Skull width----
SW.u <- Exa_data[!(Exa_data$SW.D==1),]
lm_basic_SW <- lm(SW ~ BSL, data=SW.u, offset=1.00*BSL)

I <- rep(0, length(Exa_data$SW))
I[Exa_data$SW.D==1] <- 1:sum(Exa_data$SW.D)
lmer_basic_SW <- lmer(Exa_data$SW ~ Exa_data$BSL + (1|I), offset=1.00*Exa_data$BSL)

summary(lmer_basic_SW)
summary(lm_basic_SW)

y1.sw <- (1+summary(lmer_basic_SW)$coef[2,1])*x1 + summary(lmer_basic_SW)$coef[1,1]

lmer_basic_SW <- lmer(Exa_data$SW ~ Exa_data$BSL + (1|I))
plot(Exa_data$SW ~ Exa_data$BSL, pch=16, cex=1, col="Red", 
     xlab="Skull Length (Log[mm])", ylab="Skull Width (Log[mm])")
newx.sw <- seq(min(SW.u$BSL), max(SW.u$BSL), length.out=4)
abline(lm(SW.u$SW ~ SW.u$BSL), col="black", lw=2)
conf_interval <- predict(lm_basic_SW, newdata.sw=data.frame(x=newx.sw), interval="confidence", level=0.95)
x.sw <- SW.u$BSL
lines(x.sw[order(x.sw)], conf_interval[order(x.sw),2], lty=2, col="Black", lw=2)
lines(x.sw[order(x.sw)], conf_interval[order(x.sw),3], lty=2, col="Black", lw=2)
lines(y1.sw ~ x1, col="Blue", lw=2)
m.sw <- lmer_basic_SW
newdat.sw <- data.frame(x=seq(2.2, 2.7, length.out = 23))
mm.sw <- model.matrix(~x, newdat.sw)
newdat.sw$y <- mm.sw%*%fixef(m.sw)
pvar1.sw <- diag(mm.sw %*% tcrossprod(vcov(m.sw),mm.sw))
tvar1.sw <- pvar1.sw+VarCorr(m.sw)$f[1]
newdat.sw <- data.frame(
        newdat.sw
        , plo.sw = newdat.sw$y-1.96*sqrt(pvar1.sw)
        , phi.sw = newdat.sw$y+1.96*sqrt(pvar1.sw)
)
lines(newdat.sw$x, newdat.sw$plo.sw, col="Blue", lty=2, lw=2)
lines(newdat.sw$x, newdat.sw$phi.sw, col="Blue", lty=2, lw=2)
points(SW.u$SW ~ SW.u$BSL, col="black", pch=16)

###Maxillary bicanine width----
BW.u <- Exa_data[!(Exa_data$BW.D==1),]
lm_basic_BW <- lm(BW ~ BSL, data=BW.u, offset=1.00*BSL)

I <- rep(0, length(Exa_data$BW))
I[Exa_data$BW.D==1] <- 1:sum(Exa_data$BW.D)
lmer_basic_BW <- lmer(Exa_data$BW ~ Exa_data$BSL + (1|I), offset=1.00*Exa_data$BSL)

summary(lmer_basic_BW)
summary(lm_basic_BW)

y1.bw <- (1+summary(lmer_basic_BW)$coef[2,1])*x1 + summary(lmer_basic_BW)$coef[1,1]

lmer_basic_BW <- lmer(Exa_data$BW ~ Exa_data$BSL + (1|I))
plot(Exa_data$BW ~ Exa_data$BSL, pch=16, cex=1, col="Red", 
     xlab="Skull Length (Log[mm])", ylab="Maxillary Bicanine Width (Log[mm])")
newx.bw <- seq(min(BW.u$BSL), max(BW.u$BSL), length.out=4)
abline(lm(BW.u$BW ~ BW.u$BSL), col="black", lw=2)
conf_interval <- predict(lm_basic_BW, newdata.bw=data.frame(x=newx.bw), interval="confidence", level=0.95)
x.bw <- BW.u$BSL
lines(x.bw[order(x.bw)], conf_interval[order(x.bw),2], lty=2, col="Black", lw=2)
lines(x.bw[order(x.bw)], conf_interval[order(x.bw),3], lty=2, col="Black", lw=2)
lines(y1.bw ~ x1, col="Blue", lw=2)
m.bw <- lmer_basic_BW
newdat.bw <- data.frame(x=seq(2.2, 2.7, length.out = 23))
mm.bw <- model.matrix(~x, newdat.bw)
newdat.bw$y <- mm.bw%*%fixef(m.bw)
pvar1.bw <- diag(mm.bw %*% tcrossprod(vcov(m.bw),mm.bw))
tvar1.bw <- pvar1.bw+VarCorr(m.bw)$f[1]
newdat.bw <- data.frame(
        newdat.bw
        , plo.bw = newdat.bw$y-1.96*sqrt(pvar1.bw)
        , phi.bw = newdat.bw$y+1.96*sqrt(pvar1.bw)
)
lines(newdat.bw$x, newdat.bw$plo.bw, col="Blue", lty=2, lw=2)
lines(newdat.bw$x, newdat.bw$phi.bw, col="Blue", lty=2, lw=2)
points(BW.u$BW ~ BW.u$BSL, col="black", pch=16)

###Upper postcanine tooth row----
UP.u <- Exa_data[!(Exa_data$UP.D==1),]
lm_basic_UP <- lm(UP ~ BSL, data=UP.u, offset=1.00*BSL)

I <- rep(0, length(Exa_data$UP))
I[Exa_data$UP.D==1] <- 1:sum(Exa_data$UP.D)
lmer_basic_UP <- lmer(Exa_data$UP ~ Exa_data$BSL + (1|I), offset=1.00*Exa_data$BSL)

summary(lmer_basic_UP)
summary(lm_basic_UP)

y1.up <- (1+summary(lmer_basic_UP)$coef[2,1])*x1 + summary(lmer_basic_UP)$coef[1,1]

lmer_basic_UP <- lmer(Exa_data$UP ~ Exa_data$BSL + (1|I))
plot(Exa_data$UP ~ Exa_data$BSL, pch=16, cex=1, col="Red", 
     xlab="Skull Length (Log[mm])", ylab="Upper Postcanine Tooth Row Length (Log[mm])")
newx.up <- seq(min(UP.u$BSL), max(UP.u$BSL), length.out=4)
abline(lm(UP.u$UP ~ UP.u$BSL), col="black", lw=2)
conf_interval <- predict(lm_basic_UP, newdata.up=data.frame(x=newx.up), interval="confidence", level=0.95)
x.up <- UP.u$BSL
lines(x.up[order(x.up)], conf_interval[order(x.up),2], lty=2, col="Black", lw=2)
lines(x.up[order(x.up)], conf_interval[order(x.up),3], lty=2, col="Black", lw=2)
lines(y1.up ~ x1, col="Blue", lw=2)
m.up <- lmer_basic_UP
newdat.up <- data.frame(x=seq(2.2, 2.7, length.out = 23))
mm.up <- model.matrix(~x, newdat.up)
newdat.up$y <- mm.up%*%fixef(m.up)
pvar1.up <- diag(mm.up %*% tcrossprod(vcov(m.up),mm.up))
tvar1.up <- pvar1.up+VarCorr(m.up)$f[1]
newdat.up <- data.frame(
        newdat.up
        , plo.up = newdat.up$y-1.96*sqrt(pvar1.up)
        , phi.up = newdat.up$y+1.96*sqrt(pvar1.up)
)
lines(newdat.up$x, newdat.up$plo.up, col="Blue", lty=2, lw=2)
lines(newdat.up$x, newdat.up$phi.up, col="Blue", lty=2, lw=2)
points(UP.u$UP ~ UP.u$BSL, col="black", pch=16)

###Posterior postcanine distance----
PD.u <- Exa_data[!(Exa_data$PD.D==1),]
lm_basic_PD <- lm(PD ~ BSL, data=PD.u, offset=1.00*BSL)

I <- rep(0, length(Exa_data$PD))
I[Exa_data$PD.D==1] <- 1:sum(Exa_data$PD.D)
lmer_basic_PD <- lmer(Exa_data$PD ~ Exa_data$BSL + (1|I), offset=1.00*Exa_data$BSL)

summary(lmer_basic_PD)
summary(lm_basic_PD)

y1.pd <- (1+summary(lmer_basic_PD)$coef[2,1])*x1 + summary(lmer_basic_PD)$coef[1,1]

lmer_basic_PD <- lmer(Exa_data$PD ~ Exa_data$BSL + (1|I))
plot(Exa_data$PD ~ Exa_data$BSL, pch=16, cex=1, col="Red", 
     xlab="Skull Length (Log[mm])", ylab="Posterior Postcanine Distance (Log[mm])")
newx.pd <- seq(min(PD.u$BSL), max(PD.u$BSL), length.out=4)
abline(lm(PD.u$PD ~ PD.u$BSL), col="black", lw=2)
conf_interval <- predict(lm_basic_PD, newdata.pd=data.frame(x=newx.pd), interval="confidence", level=0.95)
x.pd <- PD.u$BSL
lines(x.pd[order(x.pd)], conf_interval[order(x.pd),2], lty=2, col="Black", lw=2)
lines(x.pd[order(x.pd)], conf_interval[order(x.pd),3], lty=2, col="Black", lw=2)
lines(y1.pd ~ x1, col="Blue", lw=2)
m.pd <- lmer_basic_PD
newdat.pd <- data.frame(x=seq(2.2, 2.7, length.out = 23))
mm.pd <- model.matrix(~x, newdat.pd)
newdat.pd$y <- mm.pd%*%fixef(m.pd)
pvar1.pd <- diag(mm.pd %*% tcrossprod(vcov(m.pd),mm.pd))
tvar1.pd <- pvar1.pd+VarCorr(m.pd)$f[1]
newdat.pd <- data.frame(
        newdat.pd
        , plo.pd = newdat.pd$y-1.96*sqrt(pvar1.pd)
        , phi.pd = newdat.pd$y+1.96*sqrt(pvar1.pd)
)
lines(newdat.pd$x, newdat.pd$plo.pd, col="Blue", lty=2, lw=2)
lines(newdat.pd$x, newdat.pd$phi.pd, col="Blue", lty=2, lw=2)
points(PD.u$PD ~ PD.u$BSL, col="black", pch=16)

###Transverse process width----
TP.u <- Exa_data[!(Exa_data$TP.D==1),]
lm_basic_TP <- lm(TP ~ BSL, data=TP.u, offset=1.00*BSL)

I <- rep(0, length(Exa_data$TP))
I[Exa_data$TP.D==1] <- 1:sum(Exa_data$TP.D)
lmer_basic_TP <- lmer(Exa_data$TP ~ Exa_data$BSL + (1|I), offset=1.00*Exa_data$BSL)

summary(lmer_basic_TP)
summary(lm_basic_TP)

y1.tp <- (1+summary(lmer_basic_TP)$coef[2,1])*x1 + summary(lmer_basic_TP)$coef[1,1]

lmer_basic_TP <- lmer(Exa_data$TP ~ Exa_data$BSL + (1|I))
plot(Exa_data$TP ~ Exa_data$BSL, pch=16, cex=1, col="Red", 
     xlab="Skull Length (Log[mm])", ylab="(Transverse Process Width (Log[mm])")
newx.tp <- seq(min(TP.u$BSL), max(TP.u$BSL), length.out=4)
abline(lm(TP.u$TP ~ TP.u$BSL), col="black", lw=2)
conf_interval <- predict(lm_basic_TP, newdata.tp=data.frame(x=newx.tp), interval="confidence", level=0.95)
x.tp <- TP.u$BSL
lines(x.tp[order(x.tp)], conf_interval[order(x.tp),2], lty=2, col="Black", lw=2)
lines(x.tp[order(x.tp)], conf_interval[order(x.tp),3], lty=2, col="Black", lw=2)
lines(y1.tp ~ x1, col="Blue", lw=2)
m.tp <- lmer_basic_TP
newdat.tp <- data.frame(x=seq(2.2, 2.7, length.out = 23))
mm.tp <- model.matrix(~x, newdat.tp)
newdat.tp$y <- mm.tp%*%fixef(m.tp)
pvar1.tp <- diag(mm.tp %*% tcrossprod(vcov(m.tp),mm.tp))
tvar1.tp <- pvar1.tp+VarCorr(m.tp)$f[1]
newdat.tp <- data.frame(
        newdat.tp
        , plo.tp = newdat.tp$y-1.96*sqrt(pvar1.tp)
        , phi.tp = newdat.tp$y+1.96*sqrt(pvar1.tp)
)
lines(newdat.tp$x, newdat.tp$plo.tp, col="Blue", lty=2, lw=2)
lines(newdat.tp$x, newdat.tp$phi.tp, col="Blue", lty=2, lw=2)
points(TP.u$TP ~ TP.u$BSL, col="black", pch=16)

###Orbital diameter----
OD.u <- Exa_data[!(Exa_data$OD.D==1),]
lm_basic_OD <- lm(OD ~ BSL, data=OD.u, offset=1.00*BSL)

I <- rep(0, length(Exa_data$OD))
I[Exa_data$OD.D==1] <- 1:sum(Exa_data$OD.D)
lmer_basic_OD <- lmer(Exa_data$OD ~ Exa_data$BSL + (1|I), offset=1.00*Exa_data$BSL)

summary(lmer_basic_OD)
summary(lm_basic_OD)

y1.od <- (1+summary(lmer_basic_OD)$coef[2,1])*x1 + summary(lmer_basic_OD)$coef[1,1]

lmer_basic_OD <- lmer(Exa_data$OD ~ Exa_data$BSL + (1|I))
plot(Exa_data$OD ~ Exa_data$BSL, pch=16, cex=1, col="Red", 
     xlab="Skull Length (Log[mm])", ylab="Orbital Diameter (Log[mm])")
newx.od <- seq(min(OD.u$BSL), max(OD.u$BSL), length.out=4)
abline(lm(OD.u$OD ~ OD.u$BSL), col="black", lw=2)
conf_interval <- predict(lm_basic_OD, newdata.od=data.frame(x=newx.od), interval="confidence", level=0.95)
x.od <- OD.u$BSL
lines(x.od[order(x.od)], conf_interval[order(x.od),2], lty=2, col="Black", lw=2)
lines(x.od[order(x.od)], conf_interval[order(x.od),3], lty=2, col="Black", lw=2)
lines(y1.od ~ x1, col="Blue", lw=2)
m.od <- lmer_basic_OD
newdat.od <- data.frame(x=seq(2.2, 2.7, length.out = 23))
mm.od <- model.matrix(~x, newdat.od)
newdat.od$y <- mm.od%*%fixef(m.od)
pvar1.od <- diag(mm.od %*% tcrossprod(vcov(m.od),mm.od))
tvar1.od <- pvar1.od+VarCorr(m.od)$f[1]
newdat.od <- data.frame(
        newdat.od
        , plo.od = newdat.od$y-1.96*sqrt(pvar1.od)
        , phi.od = newdat.od$y+1.96*sqrt(pvar1.od)
)
lines(newdat.od$x, newdat.od$plo.od, col="Blue", lty=2, lw=2)
lines(newdat.od$x, newdat.od$phi.od, col="Blue", lty=2, lw=2)
points(OD.u$OD ~ OD.u$BSL, col="black", pch=16)

###Occipital plate width----
OW.u <- Exa_data[!(Exa_data$OW.D==1),]
lm_basic_OW <- lm(OW ~ BSL, data=OW.u, offset=1.00*BSL)

I <- rep(0, length(Exa_data$OW))
I[Exa_data$OW.D==1] <- 1:sum(Exa_data$OW.D)
lmer_basic_OW <- lmer(Exa_data$OW ~ Exa_data$BSL + (1|I), offset=1.00*Exa_data$BSL)

summary(lmer_basic_OW)
summary(lm_basic_OW)

y1.ow <- (1+summary(lmer_basic_OW)$coef[2,1])*x1 + summary(lmer_basic_OW)$coef[1,1]

lmer_basic_OW <- lmer(Exa_data$OW ~ Exa_data$BSL + (1|I))
plot(Exa_data$OW ~ Exa_data$BSL, pch=16, cex=1, col="Red", 
     xlab="Skull Length (Log[mm])", ylab="Occipital Plate Width (Log[mm])")
newx.ow <- seq(min(OW.u$BSL), max(OW.u$BSL), length.out=4)
abline(lm(OW.u$OW ~ OW.u$BSL), col="black", lw=2)
conf_interval <- predict(lm_basic_OW, newdata.ow=data.frame(x=newx.ow), interval="confidence", level=0.95)
x.ow <- OW.u$BSL
lines(x.ow[order(x.ow)], conf_interval[order(x.ow),2], lty=2, col="Black", lw=2)
lines(x.ow[order(x.ow)], conf_interval[order(x.ow),3], lty=2, col="Black", lw=2)
lines(y1.ow ~ x1, col="Blue", lw=2)
m.ow <- lmer_basic_OW
newdat.ow <- data.frame(x=seq(2.2, 2.7, length.out = 23))
mm.ow <- model.matrix(~x, newdat.ow)
newdat.ow$y <- mm.ow%*%fixef(m.ow)
pvar1.ow <- diag(mm.ow %*% tcrossprod(vcov(m.ow),mm.ow))
tvar1.ow <- pvar1.ow+VarCorr(m.ow)$f[1]
newdat.ow <- data.frame(
        newdat.ow
        , plo.ow = newdat.ow$y-1.96*sqrt(pvar1.ow)
        , phi.ow = newdat.ow$y+1.96*sqrt(pvar1.ow)
)
lines(newdat.ow$x, newdat.ow$plo.ow, col="Blue", lty=2, lw=2)
lines(newdat.ow$x, newdat.ow$phi.ow, col="Blue", lty=2, lw=2)
points(OW.u$OW ~ OW.u$BSL, col="black", pch=16)

###Basicranial length----
BB.u <- Exa_data[!(Exa_data$BB.D==1),]
lm_basic_BB <- lm(BB ~ BSL, data=BB.u, offset=1.00*BSL)

I <- rep(0, length(Exa_data$BB))
I[Exa_data$BB.D==1] <- 1:sum(Exa_data$BB.D)
lmer_basic_BB <- lmer(Exa_data$BB ~ Exa_data$BSL + (1|I), offset=1.00*Exa_data$BSL)

summary(lmer_basic_BB)
summary(lm_basic_BB)

y1.bb <- (1+summary(lmer_basic_BB)$coef[2,1])*x1 + summary(lmer_basic_BB)$coef[1,1]

lmer_basic_BB <- lmer(Exa_data$BB ~ Exa_data$BSL + (1|I))
plot(Exa_data$BB ~ Exa_data$BSL, pch=16, cex=1, col="Red", 
     xlab="Skull Length (Log[mm])", ylab="Basicranial Length (Log[mm])")
newx.bb <- seq(min(BB.u$BSL), max(BB.u$BSL), length.out=4)
abline(lm(BB.u$BB ~ BB.u$BSL), col="black", lw=2)
conf_interval <- predict(lm_basic_BB, newdata.bb=data.frame(x=newx.bb), interval="confidence", level=0.95)
x.bb <- BB.u$BSL
lines(x.bb[order(x.bb)], conf_interval[order(x.bb),2], lty=2, col="Black", lw=2)
lines(x.bb[order(x.bb)], conf_interval[order(x.bb),3], lty=2, col="Black", lw=2)
lines(y1.bb ~ x1, col="Blue", lw=2)
m.bb <- lmer_basic_BB
newdat.bb <- data.frame(x=seq(2.2, 2.7, length.out = 23))
mm.bb <- model.matrix(~x, newdat.bb)
newdat.bb$y <- mm.bb%*%fixef(m.bb)
pvar1.bb <- diag(mm.bb %*% tcrossprod(vcov(m.bb),mm.bb))
tvar1.bb <- pvar1.bb+VarCorr(m.bb)$f[1]
newdat.bb <- data.frame(
        newdat.bb
        , plo.bb = newdat.bb$y-1.96*sqrt(pvar1.bb)
        , phi.bb = newdat.bb$y+1.96*sqrt(pvar1.bb)
)
lines(newdat.bb$x, newdat.bb$plo.bb, col="Blue", lty=2, lw=2)
lines(newdat.bb$x, newdat.bb$phi.bb, col="Blue", lty=2, lw=2)
points(BB.u$BB ~ BB.u$BSL, col="black", pch=16)

###Zygomatic height----
ZH.u <- Exa_data[!(Exa_data$ZH.D==1),]
lm_basic_ZH <- lm(ZH ~ BSL, data=ZH.u, offset=1.00*BSL)

I <- rep(0, length(Exa_data$ZH))
I[Exa_data$ZH.D==1] <- 1:sum(Exa_data$ZH.D)
lmer_basic_ZH <- lmer(Exa_data$ZH ~ Exa_data$BSL + (1|I), offset=1.00*Exa_data$BSL)

summary(lmer_basic_ZH)
summary(lm_basic_ZH)

y1.zh <- (1+summary(lmer_basic_ZH)$coef[2,1])*x1 + summary(lmer_basic_ZH)$coef[1,1]

lmer_basic_ZH <- lmer(Exa_data$ZH ~ Exa_data$BSL + (1|I))
plot(Exa_data$ZH ~ Exa_data$BSL, pch=16, cex=1, col="Red", 
     xlab="Skull Length (Log[mm])", ylab="Zygoma Height (Log[mm])")
newx.zh <- seq(min(ZH.u$BSL), max(ZH.u$BSL), length.out=4)
abline(lm(ZH.u$ZH ~ ZH.u$BSL), col="black", lw=2)
conf_interval <- predict(lm_basic_ZH, newdata.zh=data.frame(x=newx.zh), interval="confidence", level=0.95)
x.zh <- ZH.u$BSL
lines(x.zh[order(x.zh)], conf_interval[order(x.zh),2], lty=2, col="Black", lw=2)
lines(x.zh[order(x.zh)], conf_interval[order(x.zh),3], lty=2, col="Black", lw=2)
lines(y1.zh ~ x1, col="Blue", lw=2)
m.zh <- lmer_basic_ZH
newdat.zh <- data.frame(x=seq(2.2, 2.7, length.out = 23))
mm.zh <- model.matrix(~x, newdat.zh)
newdat.zh$y <- mm.zh%*%fixef(m.zh)
pvar1.zh <- diag(mm.zh %*% tcrossprod(vcov(m.zh),mm.zh))
tvar1.zh <- pvar1.zh+VarCorr(m.zh)$f[1]
newdat.zh <- data.frame(
        newdat.zh
        , plo.zh = newdat.zh$y-1.96*sqrt(pvar1.zh)
        , phi.zh = newdat.zh$y+1.96*sqrt(pvar1.zh)
)
lines(newdat.zh$x, newdat.zh$plo.zh, col="Blue", lty=2, lw=2)
lines(newdat.zh$x, newdat.zh$phi.zh, col="Blue", lty=2, lw=2)
points(ZH.u$ZH ~ ZH.u$BSL, col="black", pch=16)

###Diastema length----
DL.u <- Exa_data[!(Exa_data$DL.D==1),]
lm_basic_DL <- lm(DL ~ BSL, data=DL.u, offset=1.00*BSL)

I <- rep(0, length(Exa_data$DL))
I[Exa_data$DL.D==1] <- 1:sum(Exa_data$DL.D)
lmer_basic_DL <- lmer(Exa_data$DL ~ Exa_data$BSL + (1|I), offset=1.00*Exa_data$BSL)

summary(lmer_basic_DL)
summary(lm_basic_DL)

y1.dl <- (1+summary(lmer_basic_DL)$coef[2,1])*x1 + summary(lmer_basic_DL)$coef[1,1]

lmer_basic_DL <- lmer(Exa_data$DL ~ Exa_data$BSL + (1|I))
plot(Exa_data$DL ~ Exa_data$BSL, pch=16, cex=1, col="Red", 
     xlab="Skull Length (Log[mm])", ylab="Diastema length (Log[mm])")
newx.dl <- seq(min(DL.u$BSL), max(DL.u$BSL), length.out=4)
abline(lm(DL.u$DL ~DL.u$BSL), col="black", lw=2)
conf_interval <- predict(lm_basic_DL, newdata.dl=data.frame(x=newx.dl), interval="confidence", level=0.95)
x.dl <- DL.u$BSL
lines(x.dl[order(x.dl)], conf_interval[order(x.dl),2], lty=2, col="Black", lw=2)
lines(x.dl[order(x.dl)], conf_interval[order(x.dl),3], lty=2, col="Black", lw=2)
lines(y1.dl ~ x1, col="Blue", lw=2)
m.dl <- lmer_basic_DL
newdat.dl <- data.frame(x=seq(2.2, 2.7, length.out=23))
mm.dl <- model.matrix(~x, newdat.dl)
newdat.dl$y <- (mm.dl%*%fixef(m.dl))
pvar1.dl <- diag(mm.dl %*% tcrossprod(vcov(m.dl),mm.dl))
tvar1.dl <- pvar1.dl+VarCorr(m.dl)$f[1]
newdat.dl <- data.frame(
        newdat.dl
        , plo.dl = newdat.dl$y-1.96*sqrt(pvar1.dl)
        , phi.dl = newdat.dl$y+1.96*sqrt(pvar1.dl)
)
lines(newdat.dl$x, newdat.dl$plo.dl, col="Blue", lty=2, lw=2)
lines(newdat.dl$x, newdat.dl$phi.dl, col="Blue", lty=2, lw=2)
points(DL.u$DL ~ DL.u$BSL, col="black", pch=16)

###Posterior zygoma width----
ZW.u <- Exa_data[!(Exa_data$ZW.D==1),]
lm_basic_ZW <- lm(ZW ~ BSL, data=ZW.u, offset=1.00*BSL)

I <- rep(0, length(Exa_data$ZW))
I[Exa_data$ZW.D==1] <- 1:sum(Exa_data$ZW.D)
lmer_basic_ZW <- lmer(Exa_data$ZW ~ Exa_data$BSL + (1|I),  offset=1.00*Exa_data$BSL)

summary(lmer_basic_ZW)
summary(lm_basic_ZW)

y1.zw <- (1+summary(lmer_basic_ZW)$coef[2,1])*x1 + summary(lmer_basic_ZW)$coef[1,1]

lmer_basic_ZW <- lmer(Exa_data$ZW ~ Exa_data$BSL + (1|I))
plot(Exa_data$ZW ~ Exa_data$BSL, pch=16, cex=1, col="Red", 
     xlab="Skull Length (Log[mm])", ylab="Zygoma Width (Log[mm])")
newx.zw <- seq(min(ZW.u$BSL), max(ZW.u$BSL), length.out=10)
abline(lm(ZW.u$ZW ~ZW.u$BSL), col="black", lw=2)
conf_interval <- predict(lm_basic_ZW, newdata.zw=data.frame(x=newx.zw), interval="confidence", level=0.95)
x.zw <- ZW.u$BSL
lines(x.zw[order(x.zw)], conf_interval[order(x.zw),2], lty=2, col="Black", lw=2)
lines(x.zw[order(x.zw)], conf_interval[order(x.zw),3], lty=2, col="Black", lw=2)
lines(y1.zw ~ x1, col="Blue", lw=2)
m.zw <- lmer_basic_ZW
newdat.zw <- data.frame(x=seq(2.2, 2.7, length.out=23))
mm.zw <- model.matrix(~x, newdat.zw)
newdat.zw$y <- (mm.zw%*%fixef(m.zw))
pvar1.zw <- diag(mm.zw %*% tcrossprod(vcov(m.zw),mm.zw))
tvar1.zw <- pvar1.zw+VarCorr(m.zw)$f[1]
newdat.zw <- data.frame(
        newdat.zw
        , plo.zw = newdat.zw$y-1.96*sqrt(pvar1.zw)
        , phi.zw = newdat.zw$y+1.96*sqrt(pvar1.zw)
)
lines(newdat.zw$x, newdat.zw$plo.zw, col="Blue", lty=2, lw=2)
lines(newdat.zw$x, newdat.zw$phi.zw, col="Blue", lty=2, lw=2)
points(ZW.u$ZW ~ ZW.u$BSL, col="black", pch=16)

#Comparing slopes for ZH and TEL
#TEL Slope = 1.1761
I <- rep(0, length(Exa_data$ZH))
I[Exa_data$ZH.D==1] <- 1:sum(Exa_data$ZH.D)
lmer_basic_ZH <- lmer(Exa_data$ZH ~ Exa_data$BSL + (1|I), offset=1.1761*Exa_data$BSL)

summary(lmer_basic_ZH)
#ZH Slope = 1.3084
I <- rep(0, length(Exa_data$TEL))
I[Exa_data$TEL.D==1] <- 1:sum(Exa_data$TEL.D)
lmer_basic_TEL <- lmer(Exa_data$TEL ~ Exa_data$BSL + (1|I), offset=1.3084*Exa_data$BSL)

summary(lmer_basic_TEL)