#=========================================================================
#=====
#===== Using log-log linear regression to model Zpol as a function of mass
#===== Citation: McGuire, Ourfalian, Ezell, Lee (in prep)
#=========================================================================
Zp_mass.df = 
  data.frame(element=rep(c("Humerus","Ulna","Radius","Femur","Tibiotarsus"),
    c(17,17,17,17,17)),
  mass=rep(c(15,83,131,242,263,209,209,341,372,314,360,364,374,455,
    498,482,563),5),
  length=c(10.9,18.0,22.4,38.3,35.4,39.9,41.1,47.4,47.8,47.5,46.7,47.6,
             49.9,49.2,50.9,50.0,51.0,
           11.4,19.5,24.3,42.8,40.1,45.4,46.2,54.1,54.5,54.5,54.1,54.9,
             56.9,56.4,56.3,58.0,56.2,
           10.0,19.4,22.5,40.3,36.1,41.1,42.4,49.0,50.3,49.5,48.3,53.4,
             50.6,55.3,52.4,54.0,52.3,
           14.0,25.4,26.5,35.2,36.8,35.7,35.4,42.5,42.4,44.1,41.7,41.5,
             44.3,44.4,45.0,45.0,45.8,
           15.2,28.7,29.8,46.6,44.2,48.4,47.6,59.5,59.9,59.3,59.7,64.3,
             61.2,60.5,62.7,61.8,63.6),
  Zp=c(0.059,0.429,0.515,4.058,2.905,4.576,5.341,10.893,12.124,11.293,
         12.818,13.649,13.123,11.510,14.294,16.442,17.822,
       0.028,0.136,0.230,2.348,1.248,2.933,3.518,5.340,6.496,5.758,6.090,
         7.029,5.649,7.016,7.058,7.797,8.053,
       0.008,0.032,0.050,0.546,0.276,0.734,0.878,1.281,1.501,1.380,1.709,
         2.158,1.481,1.845,2.205,2.266,2.516,
       0.062,0.065,0.694,2.072,2.376,2.139,2.294,4.904,4.783,4.334,5.265,
         5.162,5.282,5.544,6.196,7.558,7.668,
       0.060,0.275,0.466,1.804,1.583,1.713,1.992,3.222,4.023,3.261,4.225,
         5.311,3.545,4.095,4.759,5.681,6.444))


Zp_mass.df$size.proxy <- with(Zp_mass.df,mass)

#=====================
# Scaling of Zp
# Isometry is Zp ~ M^1
#=====================

Zp_hum.lm <- lm(log10(Zp) ~ log10(size.proxy),
  data=Zp_mass.df[Zp_mass.df$element=="Humerus",])
summary(Zp_hum.lm)
shapiro.test(Zp_hum.lm$resid) # residuals are not normal

Zp_ulna.lm <- lm(log10(Zp) ~ log10(size.proxy),
  data=Zp_mass.df[Zp_mass.df$element=="Ulna",])
summary(Zp_ulna.lm)
shapiro.test(Zp_ulna.lm$resid) # residuals are normal

Zp_rad.lm <- lm(log10(Zp) ~ log10(size.proxy),
  data=Zp_mass.df[Zp_mass.df$element=="Radius",])
summary(Zp_rad.lm)
shapiro.test(Zp_rad.lm$resid) # residuals are not normal

Zp_fem.lm <- lm(log10(Zp) ~ log10(size.proxy),
  data=Zp_mass.df[Zp_mass.df$element=="Femur",])
summary(Zp_fem.lm)
shapiro.test(Zp_fem.lm$resid) # residuals are not normal

Zp_tbt.lm <- lm(log10(Zp) ~ log10(size.proxy),
  data=Zp_mass.df[Zp_mass.df$element=="Tibiotarsus",])
summary(Zp_tbt.lm)
shapiro.test(Zp_tbt.lm$resid) # residuals are normal


#-------------------------------------------------------------------------
# Can't use t-test to test isometry b/c some residuals are not normal.
# Instead test with 95% confidence interval of scaling exponent calculated
# from non-parametric bootstrapping
#-------------------------------------------------------------------------

# Bootstrap CI function
np.boot.CI <- function(fit.res, B, beta.coef, alpha){
  fit.coef.boot <- lapply(fit.res, function(x) coef(x))
  coef.df <- t(array(unlist(fit.coef.boot), 
    dim=c(length(fit.coef.boot[[1]]),B)))
  colnames(coef.df) <- names(fit.coef.boot[[1]])
  return(quantile(coef.df[,beta.coef], probs=c(alpha/2,1-alpha/2), type=6))
}

# Bootstrap confidence bands function
np.conf.bands <- function(new.x, B, pred, alpha=0.05){
  pred.combine <- do.call(rbind, pred)
  conf.bands <- cbind(new.x, 
    t(apply(pred.combine, 2, quantile, probs=c(alpha/2, 1-alpha/2), type=6)))
  colnames(conf.bands)[2] <- "lower.band"
  colnames(conf.bands)[3] <- "upper.band"
  return(conf.bands)
}

# Resample data and repeat linear regression B times 
B <- 5000
Zp.model <- log10(Zp) ~ log10(size.proxy)
hum.orig <- Zp_mass.df[Zp_mass.df$element=="Humerus",]
ulna.orig <- Zp_mass.df[Zp_mass.df$element=="Ulna",]
rad.orig <- Zp_mass.df[Zp_mass.df$element=="Radius",]
fem.orig <- Zp_mass.df[Zp_mass.df$element=="Femur",]
tbt.orig <- Zp_mass.df[Zp_mass.df$element=="Tibiotarsus",]
pred.size <- data.frame(size.proxy=round(seq(
                 from=min(Zp_mass.df[Zp_mass.df$element=="Humerus",]$mass), 
                 to=max(Zp_mass.df[Zp_mass.df$element=="Humerus",]$mass), 
                 length.out=100),0))

hum.res <- replicate(B, hum.orig[sample(1:nrow(hum.orig),replace=T),],
  simplify=F)
hum.Zp.model.boot <- lapply(hum.res,
  function(x) lm(Zp.model, data = x))
hum.Zp.pred.res <- lapply(hum.Zp.model.boot,
  function(x) predict(x, pred.size))
hum.Zp.param.a.CI <- np.boot.CI(hum.Zp.model.boot, B, 
  beta.coef=2, alpha=0.05)
hum.Zp.param.b.CI <- np.boot.CI(hum.Zp.model.boot, B, 
  beta.coef=1, alpha=0.05)
hum.Zp.conf.bands <- 
  np.conf.bands(new.x=pred.size, B, pred=hum.Zp.pred.res)

ulna.res <- replicate(B, ulna.orig[sample(1:nrow(ulna.orig),replace=T),],
  simplify=F)
ulna.Zp.model.boot <- lapply(ulna.res,
  function(x) lm(Zp.model, data = x))
ulna.Zp.pred.res <- lapply(ulna.Zp.model.boot,
  function(x) predict(x, pred.size))
ulna.Zp.param.a.CI <- np.boot.CI(ulna.Zp.model.boot, B, 
  beta.coef=2, alpha=0.05)
ulna.Zp.param.b.CI <- np.boot.CI(ulna.Zp.model.boot, B, 
  beta.coef=1, alpha=0.05)
ulna.Zp.conf.bands <- 
  np.conf.bands(new.x=pred.size, B, pred=ulna.Zp.pred.res)

rad.res <- replicate(B, rad.orig[sample(1:nrow(rad.orig),replace=T),],
  simplify=F)
rad.Zp.model.boot <- lapply(rad.res,
  function(x) lm(Zp.model, data = x))
rad.Zp.pred.res <- lapply(rad.Zp.model.boot,
  function(x) predict(x, pred.size))
rad.Zp.param.a.CI <- np.boot.CI(rad.Zp.model.boot, B, 
  beta.coef=2, alpha=0.05)
rad.Zp.param.b.CI <- np.boot.CI(rad.Zp.model.boot, B, 
  beta.coef=1, alpha=0.05)
rad.Zp.conf.bands <- 
  np.conf.bands(new.x=pred.size, B, pred=rad.Zp.pred.res)

fem.res <- replicate(B, fem.orig[sample(1:nrow(fem.orig),replace=T),],
  simplify=F)
fem.Zp.model.boot <- lapply(fem.res,
  function(x) lm(Zp.model, data = x))
fem.Zp.pred.res <- lapply(fem.Zp.model.boot,
  function(x) predict(x, pred.size))
fem.Zp.param.a.CI <- np.boot.CI(fem.Zp.model.boot, B, 
  beta.coef=2, alpha=0.05)
fem.Zp.param.b.CI <- np.boot.CI(fem.Zp.model.boot, B, 
  beta.coef=1, alpha=0.05)
fem.Zp.conf.bands <- 
  np.conf.bands(new.x=pred.size, B, pred=fem.Zp.pred.res)

tbt.res <- replicate(B, tbt.orig[sample(1:nrow(tbt.orig),replace=T),],
  simplify=F)
tbt.Zp.model.boot <- lapply(tbt.res,
  function(x) lm(Zp.model, data = x))
tbt.Zp.pred.res <- lapply(tbt.Zp.model.boot,
  function(x) predict(x, pred.size))
tbt.Zp.param.a.CI <- np.boot.CI(tbt.Zp.model.boot, B, 
  beta.coef=2, alpha=0.05)
tbt.Zp.param.b.CI <- np.boot.CI(tbt.Zp.model.boot, B, 
  beta.coef=1, alpha=0.05)
tbt.Zp.conf.bands <- 
  np.conf.bands(new.x=pred.size, B, pred=tbt.Zp.pred.res)

#======================================================
# Collects confidence band coordinates into a single DF
#======================================================
bands.poly <- data.frame(
  element = factor(c(rep("Humerus",200), rep("Ulna",200), rep("Radius",200),
                rep("Femur",200), rep("Tibiotarsus",200))),
  bandx = c(hum.Zp.conf.bands[,1], rev(hum.Zp.conf.bands[,1]),
            ulna.Zp.conf.bands[,1], rev(ulna.Zp.conf.bands[,1]),
            rad.Zp.conf.bands[,1], rev(rad.Zp.conf.bands[,1]),
            fem.Zp.conf.bands[,1], rev(fem.Zp.conf.bands[,1]),
            tbt.Zp.conf.bands[,1], rev(tbt.Zp.conf.bands[,1])),
  bandy = 10^c(hum.Zp.conf.bands[,2], rev(hum.Zp.conf.bands[,3]),
            ulna.Zp.conf.bands[,2], rev(ulna.Zp.conf.bands[,3]),
            rad.Zp.conf.bands[,2], rev(rad.Zp.conf.bands[,3]),
            fem.Zp.conf.bands[,2], rev(fem.Zp.conf.bands[,3]),
            tbt.Zp.conf.bands[,2], rev(tbt.Zp.conf.bands[,3])))


#----------------
# Initialize plot
#----------------
library(ggplot2)
library(scales)
Zp_mass.df$element <- factor(Zp_mass.df$element, levels=c("Humerus","Ulna",
  "Radius","Femur","Tibiotarsus"), 
  labels=c("Humerus","Ulna","Radius","Femur","Tibiotarsus"))
dev.new(width=4.25, height=11.2)

ggplot(Zp_mass.df, aes(x=size.proxy,y=Zp)) + 
  facet_wrap(~element, ncol=1) +
  theme_bw() + 
  theme(plot.background=element_blank(),
    panel.grid.major=element_blank(),
    panel.grid.minor=element_blank(),
    panel.border=element_rect(color="black", size=1),
    panel.background=element_blank(),
    axis.line.x=element_line(size=0),
    axis.line.y=element_line(size=0),
    axis.ticks=element_line(size=1),
    axis.ticks.length=unit(2,"mm"),
    strip.background = element_blank(),
    strip.text.x = element_blank(),
    legend.position="none") + 
  xlab(expression(paste("log"[10]," Body mass (g)"))) + 
  ylab(expression(paste("log"[10]," Polar section modulus ( ",mm^{3},")"))) + 
  scale_x_log10(breaks = trans_breaks("log10",function(x) 10^x),
    labels = trans_format("log10", math_format(.x))) + 
  scale_y_log10(breaks = trans_breaks("log10",function(x) 10^x),
    labels = trans_format("log10", math_format(.x))) +
  stat_smooth(method=lm,aes(fill=factor(element),color=factor(element)),
    size=1.25, se=FALSE, fullrange=TRUE, show.legend=FALSE) + 
  scale_color_manual(values = c("green4", "deepskyblue", "red3", "purple3",
    "orange4"), guide=guide_legend(title=NULL)) + 
  geom_point(aes(size=factor(element), fill=factor(element), 
    shape=factor(element))) + 
  scale_size_manual(values = c(3,4,4.5,4,3), guide=guide_legend(title=NULL)) +
  scale_fill_manual(values = c("green", "deepskyblue", "red", "purple", 
    "orange"), guide=guide_legend(title=NULL)) +
  scale_shape_manual(values = c(24, 22, 21, 23, 25), 
    guide=guide_legend(title=NULL)) + 
  geom_abline(intercept=-5, slope=1, size=1) +
  annotate("text", x = 360, y = 0.001, label="Isometry", size=4, angle=10) +
  annotate("text", x = 16, y = 18, label="A", size=6, fontface="bold") +
  geom_polygon(data=bands.poly, aes(x=bandx, y=bandy, group=element,
    fill=factor(element), alpha=0.5))