###Code used by Braun, Wilkinson, and Cade (2020. Comment on "Female ###toads engaging in adaptive hybridization prefer high-quality ###heterospecifics as mates”, Science.  A technical comment on Chen and ###Pfennig 2020 analyses of tadpole hybidization.

###Data set DataS1 available for download from Supplementary Material https://science.sciencemag.org/content/suppl/2020/03/18/367.6484.1377.DC1

###Code written by B.S. Cade, USGS (cadeb@usgs.gov).

tadpole.hof <- read.csv("c:\\cadestuff\\data\\MBraun\\aaz5109_DataS1.csv",header=T)

###Some weird problem with name on first column
colnames(tadpole.hof)[1] <- "BreedingFamily"

###Appears that they use prcomp() with scale = F (default)
###Note that we have sign reversal on some coefficients in PC2 and PC3 ###compared to their supplementary table output.  However, proportions ###of variance are identical so these PC1 scores if not identical to C&P ###are very close.

pca.model <- prcomp(~TadMass + TadSVL + TadGosnerStage,data=tadpole.hof)

tadpole.hof$PC1 <- pca.model$x[,1]

###C&P method
tadpole.hof$PC1.alt <- predict(pca.model,newdata=tadpole.hof)[,1]

###Used identical values for PC1.

summary(tadpole.hof$PC1 - tadpole.hof$PC1.alt)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      0       0       0       0       0       0 

###Using Philip Dixon suggestion of PCA on correlations because of huge ###differences in scales of variables. This results in PC1 only ###explaining 88% of variance.

pca.cor.model <- princomp(~TadMass + TadSVL + TadGosnerStage,data=tadpole.hof,cor=T)

summary(pca.cor.model)

tadpole.hof$PC1.cor <- pca.cor.model$scores[,1]

table(tadpole.hof$BreedingFamily)

###BreedingFamily 11 only had 7 tadpoles and 15 only had 1 tadpole

###No need to standardize or center predictors for interpreting partial ###effects of a given predictor so did not do that in our analyses.
###Model averaging of regression coefficients also not used as it is ###wrong without standardizing by partial standard deviations and is not ###really useful even then (see Cade 2015).

library(lme4)

mixed.1 <- lmer(PC1 ~ CallPulseRate + FemaleSVL + MaleCondition + (1|BreedingFamily), data=tadpole.hof)

summary(mixed.1)

mixed.2 <- lmer(PC1 ~ CallPulseRate + FemaleSVL + CallDuration + (1|BreedingFamily), data=tadpole.hof)

summary(mixed.2)

mixed.3 <- lmer(PC1 ~ CallPulseRate + FemaleSVL + (1|BreedingFamily), data=tadpole.hof)

summary(mixed.3)

mixed.4 <- lmer(PC1 ~ CallPulseRate + MaleCondition + (1|BreedingFamily), data=tadpole.hof)

summary(mixed.4)

mixed.5 <- lmer(PC1 ~ CallPulseRate + (1|BreedingFamily), data=tadpole.hof)

summary(mixed.5)

library(MuMIn)

AICc(mixed.1,mixed.2,mixed.3,mixed.4,mixed.5)

###Note that the ordering of AICc here differ from those of C&P because ###the log likelihoods calculated with the standardized predictors used ###by C&P differed from mine except for mixed.5.  There is no conceptual ###reason that this should be the case, so likely has to do with the ###numerics of mixed-effects estimation of variances with standardized ###versus unstandardized data.

###Betas for CallPulseRate for the five models
> fixef(mixed.1)[2]
CallPulseRate 
    -1.072143 
> fixef(mixed.2)[2]
CallPulseRate 
    -1.213705 
> fixef(mixed.3)[2]
CallPulseRate 
    -1.021186 
> fixef(mixed.4)[2] 
CallPulseRate 
   -0.9585506 
> fixef(mixed.5)[2] 
CallPulseRate 
   -0.9278959 
>

library(car)
###Notice how increasing VIF with increasingly negative estimate of ###CallPulseRate from -0.928 when no collinearity.

> vif(mixed.1)[1]
CallPulseRate 
     1.020832 
> vif(mixed.2)[1]
CallPulseRate 
     1.117795 
> vif(mixed.3)[1]
CallPulseRate 
     1.014178 
> vif(mixed.4)[1]
CallPulseRate 
     1.002484 
> 

###Coefficient of determinations
library(MuMIn)
r.squaredGLMM(mixed.1)
r.squaredGLMM(mixed.2)
r.squaredGLMM(mixed.3)
r.squaredGLMM(mixed.4)
r.squaredGLMM(mixed.5)

###The estimates below excluding call pulse rate can be used to get ####partial R2m by subtraction.

mixed.1.exc.pulse <- lmer(PC1 ~ FemaleSVL + MaleCondition + (1|BreedingFamily), data=tadpole.hof)
mixed.2.exc.pulse <- lmer(PC1 ~ FemaleSVL + CallDuration + (1|BreedingFamily), data=tadpole.hof)
mixed.3.exc.pulse <- lmer(PC1 ~ FemaleSVL + (1|BreedingFamily), data=tadpole.hof)
mixed.4.exc.pulse <- lmer(PC1 ~ MaleCondition + (1|BreedingFamily), data=tadpole.hof)

r.squaredGLMM(mixed.1)[1] - r.squaredGLMM(mixed.1.exc.pulse)[1]
r.squaredGLMM(mixed.2)[1]- r.squaredGLMM(mixed.2.exc.pulse)[1]
r.squaredGLMM(mixed.3)[1]- r.squaredGLMM(mixed.3.exc.pulse)[1]
r.squaredGLMM(mixed.4)[1]- r.squaredGLMM(mixed.4.exc.pulse)[1]

###R2m for fixed effects vary from 0.048 (mixed.5) to 0.119 (mixed.1).
###R2c for fixed + random effects vary from 0.258 (mixed.1) to 0.278 ###(mixed.5).  But note that R2c for model with just random effects for ###family is 0.289 with estimate of sigma = 2.528.

###Also note that since there is a single unique value of call pulse ###rate (as well as for other parental predictors) for each family, that ###treating family as a categorical variable with random effects is ###exactly the same as treating call pulse rate (other predictors) as a ###categorical predictor.  So the proportion of variance explained by ###family random effect (general functional form) compared to proportion ###explained by linear predictor of call pulse rate can be used to ###indicate lack-of-fit of the linear relationship.

mixed.family <- lmer(PC1 ~ (1|BreedingFamily), data=tadpole.hof)

summary(mixed.family)
r.squaredGLMM(mixed.family)

mixed.callpulserate <- lmer(PC1 ~ (1|CallPulseRate), data=tadpole.hof)

summary(mixed.callpulserate)
r.squaredGLMM(mixed.callpulserate)

###This yields the estimated means by family (or call pulse rate)

fixef(mixed.family) + ranef(mixed.family)$BreedingFamily

cat.means <- data.frame(fixef(mixed.callpulserate) + ranef(mixed.callpulserate)$CallPulseRate)

colnames(cat.means)[1] <- "Mean.PC1"

cat.means <- data.frame(cat.means$Mean.PC1,CallPulseRate=c(23.94287305,24.46950577,25.49138056,25.53675549,25.66271065,25.66475039,25.6895461,25.74036123,25.86866823,25.98971231,26.01044951,26.11098611,26.14152267,26.77671045,26.84483094,26.89168001,28.25336673))

colnames(cat.means)[1] <- "Mean.PC1"

###Note that including call pulse rate and other parental predictors as ###linear functions (fixed effects) barely reduces the estimated ###standard deviation of the regression (sigma = 2.526 - 2.529) compared ###to treating them as categorical predictors with Breeding Family ###(2.528). This indicates linear function of predictors not really very ###explanatory.

sigma(mixed.family)
sigma(mixed.1)
sigma(mixed.2)
sigma(mixed.3)
sigma(mixed.4)
sigma(mixed.5)

###Graphing lack-of-fit between linear function of call pulse rate and ###categorical (general function) means.

###family 15 has ranef = 0.0126 so close to zero 
x.graph <- data.frame(CallPulseRate = seq(min(tadpole.hof$CallPulseRate), max(tadpole.hof$CallPulseRate),by=0.1),BreedingFamily = c(15))

plot(PC1~CallPulseRate,data=tadpole.hof,pch=16,col="gray",xlab="Sire  pulse rate",ylab="Tadpole PC1 score",cex.lab=1.5,ylim=c(-14,6),xlim=c(23.8,28.4))
lines(cat.means$CallPulseRate,cat.means$Mean.PC1,col="red",lwd=2)
lines(x.graph$CallPulseRate,predict(mixed.5, newdata=x.graph),col="black",lwd=2)

###Using correlation structure in PCA reduces spread of PC1
plot(PC1~CallPulseRate,data=tadpole.hof)
points(PC1.cor~CallPulseRate,data=tadpole.hof,col="red")

###Alternative method examining lack-of-fit using b-splines, piecewise ###linear.
library(splines)

lm.call.bs <- lm(PC1~bs(CallPulseRate,degree=1,df=8),data=tadpole.hof)
x.graph <- data.frame(CallPulseRate = seq(min(tadpole.hof$CallPulseRate), max(tadpole.hof$CallPulseRate),by=0.1))

plot(PC1~CallPulseRate,data=tadpole.hof)
lines(x.graph$CallPulseRate,predict(lm.call.bs,newdata=x.graph),col="red")

mixed.5.bs <- lmer(PC1 ~ bs(CallPulseRate,degree=1,df=9) + (1|BreedingFamily), data=tadpole.hof)

###bs(degree=1,df=9) provides 8 knots and 2 boundary values breaking up ###CallPulserRate into 11.1%tiles

summary(mixed.5.bs)
AICc(mixed.5.bs)
AICc(mixed.5)

AICc(mixed.5) - AICc(mixed.5.bs)
[1] 14.51199
###Piecewise linear model is better fit

###family 15 has ranef = 0.0126 so close to zero 
x.graph <- data.frame(CallPulseRate = seq(min(tadpole.hof$CallPulseRate), max(tadpole.hof$CallPulseRate),by=0.1),BreedingFamily = c(15))

plot(PC1~CallPulseRate,data=tadpole.hof,pch=16,col="gray",xlab="Sire pulse rate",ylab="Tadpole PC1 score",cex.lab=1.5)
lines(x.graph$CallPulseRate,predict(mixed.5.bs,newdata=x.graph),col="red",lwd=2)
lines(x.graph$CallPulseRate,predict(mixed.5,newdata=x.graph),col="black",lwd=2)

###Deleting family 9 which seems to be driving relationship, but only ###for mixed.3 to mixed.5

mixed.1.del.fam9 <- lmer(PC1 ~ CallPulseRate + FemaleSVL + MaleCondition + (1|BreedingFamily), data=tadpole.hof[tadpole.hof$BreedingFamily!=9,])

summary(mixed.1.del.fam9)

mixed.2.del.fam9 <- lmer(PC1 ~ CallPulseRate + FemaleSVL + CallDuration + (1|BreedingFamily), data=tadpole.hof[tadpole.hof$BreedingFamily!=9,])

summary(mixed.2.del.fam9)

mixed.3.del.fam9 <- lmer(PC1 ~ CallPulseRate + FemaleSVL + (1|BreedingFamily), data=tadpole.hof[tadpole.hof$BreedingFamily!=9,])

summary(mixed.3.del.fam9)

mixed.4.del.fam9 <- lmer(PC1 ~ CallPulseRate + MaleCondition + (1|BreedingFamily), data=tadpole.hof[tadpole.hof$BreedingFamily!=9,])

summary(mixed.4.del.fam9)

mixed.5.del.fam9 <- lmer(PC1 ~ CallPulseRate + (1|BreedingFamily), data=tadpole.hof[tadpole.hof$BreedingFamily!=9,])

summary(mixed.5.del.fam9)

###
library(car)
vif(mixed.1)
vif(mixed.1.del.fam9)
###deleting family 9 increases multicollinearity between call pulse ###rate, female SVL, and Male Condition.

(vif(mixed.1.del.fam9))^-0.5
CallPulseRate     FemaleSVL MaleCondition 
    0.8833072     0.9149745     0.9055804 
vif(mixed.1))^-0.5
CallPulseRate     FemaleSVL MaleCondition 
    0.9897439     0.9822548     0.9895229 

vif(mixed.2)
vif(mixed.2.del.fam9)

plot(CallPulseRate~FemaleSVL,data=tadpole.hof)
points(CallPulseRate~FemaleSVL,data=tadpole.hof[tadpole.hof$BreedingFamily!=9,],pch=16)

plot(CallPulseRate~MaleCondition,data=tadpole.hof)
points(CallPulseRate~MaleCondition,data=tadpole.hof[tadpole.hof$BreedingFamily!=9,],pch=16)

plot(FemaleSVL~MaleCondition,data=tadpole.hof)
points(FemaleSVL~MaleCondition,data=tadpole.hof[tadpole.hof$BreedingFamily!=9,],pch=16)

plot(CallPulseRate~CallDuration,data=tadpole.hof)
points(CallPulseRate~CallDuration,data=tadpole.hof[tadpole.hof$BreedingFamily!=9,],pch=16)

###Same lack-of-fit with b-splines but for tadpole mass
mixed.5.mass.bs <- lmer(TadMass ~ bs(CallPulseRate,degree=1,df=9) + (1|BreedingFamily), data=tadpole.hof)

plot(TadMass~CallPulseRate,data=tadpole.hof,pch=16,col="gray",xlab="Sire call pulse rate",ylab="Tadpole mass(g)",main="Piecewise linear (red) and linear (black) models",cex.lab=1.5,cex.main=1.5)
lines(x.graph$CallPulseRate,predict(mixed.5.mass.bs,newdata=x.graph),col="red",lwd=2)
lines(x.graph$CallPulseRate,predict(mixed.5.mass,newdata=x.graph),col="black",lwd=2)

###This term from Bring (1994) and Cade (2015) show how much the ###residualizing effect of multicollinearity has had on the predictor ###CallPulseRate.  We will use these to alter the distance between min ###and max CallPulseRate in separation of distributions calculations ###(Xie and Nelson 2003).
 
residualize.1 <- (vif(mixed.1)[1])^-0.5
residualize.2 <- (vif(mixed.2)[1])^-0.5
residualize.3 <- (vif(mixed.3)[1])^-0.5
residualize.4 <- (vif(mixed.4)[1])^-0.5

###Separation of distributions Xie and Nelson 2003

max(tadpole.hof$CallPulseRate) - min(tadpole.hof$CallPulseRate)
[1] 4.310494

###min to max pulse rate with no collinearity (mixed.5) is 4.31 

###Quantile split QS for 4.31 between min and max CallPulseRate;
###normal distribution function((-|estimated slope for CallPulseRate| * ###separation distance)/(2*estimated SD for regression))
###Here separation distance between min and max is reduced slighlty from ###4.31 when there is multicollinearity in the residualizing effect.

pnorm(-1*(abs(fixef(mixed.1)[2]) * (4.31*residualize.1))/(2*sigma(mixed.1)))
pnorm(-1*(abs(fixef(mixed.2)[2])*(4.31*residualize.2))/(2*sigma(mixed.2)))
pnorm(-1*(abs(fixef(mixed.3)[2])*(4.31*residualize.3))/(2*sigma(mixed.3)))
pnorm(-1*(abs(fixef(mixed.4)[2])*(4.31*residualize.4))/(2*sigma(mixed.4)))
pnorm(-1*(abs(fixef(mixed.5)[2])*4.31)/(2*sigma(mixed.5)))

###Quantile splits vary from 0.16 to 0.21

CallPulseRate 
    0.1827738 
CallPulseRate 
    0.1636688
CallPulseRate 
    0.1935113 
CallPulseRate 
    0.2072692 
CallPulseRate 
    0.2144369 

###Quantile split QS for 1.93 between min and median CallPulseRate;
###normal distribution function((-|estimated slope for CallPulseRate| * ###separation distance)/(2*estimated SD for regression))
###Here separation distance between min and max is reduced slighlty from ###1.93 when there is multicollinearity in the residualizing effect.

pnorm(-1*(abs(fixef(mixed.1)[2]) * (1.93*residualize.1))/(2*sigma(mixed.1)))
pnorm(-1*(abs(fixef(mixed.2)[2])*(1.93*residualize.2))/(2*sigma(mixed.2)))
pnorm(-1*(abs(fixef(mixed.3)[2])*(1.93*residualize.3))/(2*sigma(mixed.3)))
pnorm(-1*(abs(fixef(mixed.4)[2])*(1.93*residualize.4))/(2*sigma(mixed.4)))
pnorm(-1*(abs(fixef(mixed.5)[2])*1.93)/(2*sigma(mixed.5)))

###Quantile splits vary from 0.33 to 0.36

CallPulseRate 
    0.3426705 
CallPulseRate 
    0.3304714 
CallPulseRate 
    0.3492461 
CallPulseRate 
    0.3574171 
CallPulseRate 
    0.3615719 

###Similar quantile split QS for 1.5 for majority of CallPulseRates

pnorm(-1*(abs(fixef(mixed.1)[2])*(1.5*residualize.1))/(2*sigma(mixed.1)))
pnorm(-1*(abs(fixef(mixed.2)[2])*(1.5*residualize.2))/(2*sigma(mixed.2)))
pnorm(-1*(abs(fixef(mixed.3)[2])*(1.5*residualize.3))/(2*sigma(mixed.3)))
pnorm(-1*(abs(fixef(mixed.4)[2])*(1.5*residualize.4))/(2*sigma(mixed.4)))
pnorm(-1*(abs(fixef(mixed.5)[2])*1.5)/(2*sigma(mixed.5)))

###Quantile splits vary from 0.37 to 0.39

CallPulseRate 
    0.3764145 
CallPulseRate 
    0.3665932 
CallPulseRate 
    0.3816865 
CallPulseRate 
    0.3882177 
CallPulseRate 
    0.3915306

###Graphing quantile splits for model mixed.5 with no collinearity
plot(PC1 ~ CallPulseRate,data=tadpole.hof,ylim=c(-14,6),xlim=c(23.8,28.4),ylab="Tadpole PC1 score",xlab="Sire pulse rate",cex.lab=1.5,pch=16,col="gray")
x.graph <- seq(23.8,28.4,by=0.02)

lines(x.graph,fixef(mixed.5)[1] + (fixef(mixed.5)[2])*x.graph,col="black",lwd=2)

lines(x.graph,(fixef(mixed.5)[1] + qnorm(0.361)*(sigma(mixed.5))) + (fixef(mixed.5)[2])*x.graph,col="red",lwd=1)

lines(x.graph,(fixef(mixed.5)[1] + qnorm(0.639)*(sigma(mixed.5))) + (fixef(mixed.5)[2])*x.graph,col="red",lwd=1)

lines(x.graph,(fixef(mixed.5)[1] + qnorm(0.214)*(sigma(mixed.5))) + (fixef(mixed.5)[2])*x.graph,col="blue",lwd=1)

lines(x.graph,(fixef(mixed.5)[1] + qnorm(0.796)*(sigma(mixed.5))) + (fixef(mixed.5)[2])*x.graph,col="blue",lwd=1)

###To draw flat separation of distance lines

low.361 <- (fixef(mixed.5)[1] + qnorm(0.361)*(sigma(mixed.5))) + (fixef(mixed.5)[2])*(min(tadpole.hof$CallPulseRate))

x1.temp <- c(min(tadpole.hof$CallPulseRate), min(tadpole.hof$CallPulseRate)+1.93)
y1.temp <- rep(low.361,2)

low.214 <- (fixef(mixed.5)[1] + qnorm(0.214)*(sigma(mixed.5))) + (fixef(mixed.5)[2])*(min(tadpole.hof$CallPulseRate))

x2.temp <- c(min(tadpole.hof$CallPulseRate), min(tadpole.hof$CallPulseRate)+4.31)
y2.temp <- rep(low.214,2)

lines(x1.temp,y1.temp,lwd=2,lty=2,col="red")
lines(x2.temp,y2.temp,lwd=2,lty=2,col="blue")

###Tadpole Survival

###Binomial confidence intervals on survival
tadpole.hof$numTadpoles[tadpole.hof$family==1] <- 34 
tadpole.hof$numTadpoles[tadpole.hof$family==2] <- 37 
tadpole.hof$numTadpoles[tadpole.hof$family==3] <- 35  
tadpole.hof$numTadpoles[tadpole.hof$family==4] <- 13 
tadpole.hof$numTadpoles[tadpole.hof$family==5] <- 43
tadpole.hof$numTadpoles[tadpole.hof$family==6] <- 40 
tadpole.hof$numTadpoles[tadpole.hof$family==7] <- 18 
tadpole.hof$numTadpoles[tadpole.hof$family==9] <- 15  
tadpole.hof$numTadpoles[tadpole.hof$family==11] <- 7 
tadpole.hof$numTadpoles[tadpole.hof$family==12]<- 41 
tadpole.hof$numTadpoles[tadpole.hof$family==13] <- 23 
tadpole.hof$numTadpoles[tadpole.hof$family==15] <- 1 
tadpole.hof$numTadpoles[tadpole.hof$family==16] <- 48 
tadpole.hof$numTadpoles[tadpole.hof$family==17] <- 29 
tadpole.hof$numTadpoles[tadpole.hof$family==18] <- 9 
tadpole.hof$numTadpoles[tadpole.hof$family==19] <- 43 
tadpole.hof$numTadpoles[tadpole.hof$family==20] <- 43

library(Hmisc)
###Here adding in 0 survival for three families (8,10,14) not reported
temp <- data.frame(rbind(binconf(34,50,method="exact",alpha=0.05),binconf(37,50,method="exact",alpha=0.05),binconf(35,50,method="exact",alpha=0.05),binconf(13,50,method="exact",alpha=0.05),binconf(43,50,method="exact",alpha=0.05),binconf(40,50,method="exact",alpha=0.05),binconf(18,29,method="exact",alpha=0.05),binconf(15,50,method="exact",alpha=0.05),binconf(7,9,method="exact",alpha=0.05),binconf(41,50,method="exact",alpha=0.05),binconf(23,37,method="exact",alpha=0.05),binconf(1,4,method="exact",alpha=0.05),binconf(48,50,method="exact",alpha=0.05),binconf(29,50,method="exact",alpha=0.05),binconf(9,50,method="exact",alpha=0.05),binconf(43,50,method="exact",alpha=0.05),binconf(43,50,method="exact",alpha=0.05),binconf(0,50,method="exact",alpha=0.05),binconf(0,50,method="exact",alpha=0.05),binconf(0,50,method="exact",alpha=0.05)))
summary.qr.hof <- cbind(summary.qr.hof,temp[-12,])
###Data from Table S2, including family = 15 with only a single tadpole.
SirePulseRate <- c(25.74,25.66,26.89,23.94,25.99,26.14,26.11,28.25,25.49,25.66,25.54,24.47,26.84,26.78,26.01,25.87,25.69,24.12,25.74,23.43)
TadpoleSurvival <- c(0.68,0.74,0.70,0.26,0.86,0.80,0.62,0.30,0.78,0.82,0.62,0.25,0.96,0.58,0.18,0.86,0.86,0.0,0.0,0.0)
Family <- c(1,2,3,4,5,6,7,9,11,12,13,15,16,17,18,19,20,8,10,14)
summary.hof.plus15 <- data.frame(cbind(Family,SirePulseRate,TadpoleSurvival))
summary.hof.plus15 <- cbind(summary.hof.plus15,temp)
###There is a weak, statistically nonzero, quadratic relationship
###between tadpole survival rate and sire pulse rate, driven by the ###three low survival rates at the lowest and highest pulse rates.

library(splines)

plot(TadpoleSurvival ~ SirePulseRate,data=summary.hof.plus15,ylim=c(0,1),xlim=c(23.5,28.5),ylab="Survival rate",xlab="Sire pulse rate",cex.lab=1.5,pch=16,cex=1.5)
x.graph <- data.frame(SirePulseRate = seq(min(summary.qr.hof$SirePulseRate), max(summary.qr.hof$SirePulseRate),by=0.1))

lines(x.graph$SirePulseRate,predict(tad.surv.bs.plus15,newdata=x.graph),col="black",lwd=2,ylim=c(0.0,1.0))
points(summary.hof.plus15$SirePulseRate[summary.hof.plus15$TadpoleSurvival==0.0],summary.hof.plus15$PointEst[summary.hof.plus15$TadpoleSurvival==0.0],pch=16,col="red",cex=1.5)

###Estimates with binomial confidence intervals

x.graph <- data.frame(SirePulseRate = seq(min(summary.qr.hof$SirePulseRate), max(summary.qr.hof$SirePulseRate),by=0.1))

xYplot(Cbind(PointEst,Lower,Upper)~SirePulseRate, data=summary.hof.plus15[summary.hof.plus15$TadpoleSurvival>0,],col=c("black"),ylim=c(0.0,1.0),xlim=c(23.5,28.5),pch=16,ylab=list(label="Survival rate",cex=1.5),xlab=list(label="Sire pulse rate",cex=1.5),main=list(label="Tadpole survival rates and 95% CI",cex=1.5))
###Polynomial of survival on sire pulse rate, family 15 with single ###tadpole included but families with zero survival excluded.

tad.surv.bs.plus15 <-lm(TadpoleSurvival ~ bs(SirePulseRate,degree=2,df=2),data=summary.hof.plus15[summary.hof.plus15$TadpoleSurvival>0.0,])
summary(tad.surv.bs.plus15)
###P = 0.0142 for joint effect of b-spline terms
###Drop in AICc from linear model to b-spline model = 6.35

AICc(tad.surv.bs.plus15)

AICc(lm(TadpoleSurvival ~ SirePulseRate,data=summary.hof.plus15[summary.hof.plus15$TadpoleSurvival>0.0,]))

###Polynomial of survival on sire pulse rate, family 15 with single ###tadpole deleted.
plot(TadpoleSurvival~SirePulseRate, data=summary.qr.hof,col=c("black"),ylim=c(0.0,1.0),pch=16,ylab="Tadpole survival rate",xlab="Sire pulse rate",cex.lab=1.5,main="Tadpole Survival",cex.main=1.5)
library(Hmisc)
xYplot(Cbind(PointEst,Lower,Upper)~SirePulseRate, data=summary.qr.hof,col=c("black"),ylim=c(0.0,1.0),pch=16,ylab=list(label="Survival rate",cex=1.5),xlab=list(label="Sire pulse rate",cex=1.5),main=list(label="Tadpole survival rates and 95% CI",cex=1.5))
###Excluding family 15 with single tadpole
library(splines)
summary(lm(TadpoleSurvival ~ bs(SirePulseRate,degree=2,df=2),data=summary.hof))
library(MuMIn)
summary(lm(TadpoleSurvival ~ SirePulseRate,data=summary.hof))
tad.surv.bs <-lm(TadpoleSurvival ~ bs(SirePulseRate,degree=2,df=2),data=summary.hof)
AICc(lm(TadpoleSurvival ~ bs(SirePulseRate,degree=2,df=2),data=summary.hof))
AICc(lm(TadpoleSurvival ~ bs(SirePulseRate,degree=2,df=2),data=summary.hof))
###Alternative model treating tadpole mass at SVL as measure of ###condition rather than artifical PC1 scores construct. See Cade et al. ###(2011) and Lleonart et al.(2000).

plot(TadMass~TadSVL,data=tadpole.hof,col=BreedingFamily)

###looks like typical allometric M ~ B0 SVL^B1

###Here plotting tadpole mass versus SVL for 16 families (family 15 with ###single tadpole excluded) and showing estimated median allometric ###relationships for families 4 and 9 with lowest and highest ###CallPulseRate. 

###0.50 (median) quantile regression allowing separate intercepts and ###slopes for each of the 16 families (fixed effect).

tadpole.hof$family <- tadpole.hof$BreedingFamily
tadpole.hof$family <- as.factor(tadpole.hof$family)

library(quantreg)
qr.family.het <- rq(log10(TadMass) ~ log10(TadSVL) + family + family:log10(TadSVL),tau=c(0.50),data=tadpole.hof[tadpole.hof$BreedingFamily!=15,],contrasts=list(family="contr.treatment"))

summary(qr.family.het,se="boot",bsmethod="wxy",R=1000)
###Model with just intercept
qr.family.null <- rq(log10(TadMass) ~ 1,tau=c(0.50), data=tadpole.hof[tadpole.hof$BreedingFamily!=15,],contrasts=list(family="contr.treatment"))

###Quantile regression coefficient of determination
1 - (qr.family.het$rho/qr.family.null$rho)
[1] 0.6894558
###The simpler model with a common slope (allometric exponent) for each ###family.

qr.family.common.slope <- rq(log10(TadMass) ~ log10(TadSVL) + family,tau=c(0.50),data=tadpole.hof[tadpole.hof$BreedingFamily!=15,],contrasts=list(family="contr.treatment"))

summary(qr.family.common.slope,se="boot",bsmethod="wxy",R=1000)
AICc(qr.family.het) - AICc(qr.family.common.slope)
###Difference in AICc of -32 indicates more complex model with separate ###slopes for each family is supported.
anova(qr.family.het,qr.family.common.slope,test="Wald",se="boot",iid=F,R=1000)

###Hypothesis test also supports more complex model.

###Adding in tadpole gosner stage did not improve model much (P = 0.079)

qr.family.het.gosner <- rq(log10(TadMass) ~ log10(TadSVL) + family + family:log10(TadSVL) + log10(TadGosnerStage),tau=c(0.50),data=tadpole.hof[tadpole.hof$BreedingFamily!=15,],contrasts=list(family="contr.treatment"))

summary(qr.family.het.gosner,se="boot",bsmethod="wxy",R=1000)
1 - (qr.family.het.gosner$rho/qr.family.null$rho)
[1] 0.6914345
###So including Gosner stage with P = 0.08 only increased coefficient of ###determination (variation explained) by 0.002. 
###Graphs of allometric relationships

plot(TadMass~TadSVL,data=tadpole.hof[tadpole.hof$family!=4 & tadpole.hof$family!=9 & tadpole.hof$family!=15,],col="gray",pch=16,ylab="Tadpole mass (g)",xlab="Tadpole SVL (mm)",cex.lab=1.5)
points(TadMass~TadSVL,data=tadpole.hof[tadpole.hof$family==4,],col="red",pch=16)
points(TadMass~TadSVL,data=tadpole.hof[tadpole.hof$family==9,],col="blue",pch=16)

xTadSVL <- seq(from=1,to = 17.0,by=0.1) 
lines(xTadSVL, 10^(qr.family.het$coef[1] + qr.family.het$coef[5] + (qr.family.het$coef[2] + qr.family.het$coef[20])*log10(xTadSVL)),col="red",lwd=2)
lines(xTadSVL, 10^(qr.family.het$coef[1] + qr.family.het$coef[9] + (qr.family.het$coef[2] + qr.family.het$coef[24])*log10(xTadSVL)),col="blue",lwd=2)

###Similar model estimated with mean regression. Backtransformed ###estimates are geometric mean estimates of median regression. Not used ###by us because more sensitive to a few outlying values. Also P-values ###not adjusted for heterogeneous variances.

ols.family.het <- lm(log10(TadMass) ~ log10(TadSVL) + family + family:log10(TadSVL), data=tadpole.hof[tadpole.hof$BreedingFamily!=15,],contrasts=list(family="contr.treatment"))

summary(ols.family.het)
###Getting predicted mass at specific SVL and 95% CI by families for ###comparisons. Note SVL = 14 mm is median SVL across all tadpoles in ###all families.

PredMass.SVL.14.qr50 <- 10^(predict(qr.family.het, newdata=data.frame(family=factor(c(1,2,3,4,5,6,7,9,11,12,13,16,17,18,19,20)),TadSVL=14),interval="confidence",level=0.95,type="percentile",se="boot",bsmethod="wxy",R=1000))
PredMass.SVL.12.qr50 <- 10^(predict(qr.family.het, newdata=data.frame(family=factor(c(1,2,3,4,5,6,7,9,11,12,13,16,17,18,19,20)),TadSVL=12),interval="confidence",level=0.95,type="percentile",se="boot",bsmethod="wxy",R=1000))
PredMass.SVL.16.qr50 <- 10^(predict(qr.family.het, newdata=data.frame(family=factor(c(1,2,3,4,5,6,7,9,11,12,13,16,17,18,19,20)),TadSVL=16),interval="confidence",level=0.95,type="percentile",se="boot",bsmethod="wxy",R=1000))
PredMass.SVL.14.qr50.wtd <- 10^(predict(qr.family.het.wtd, newdata=data.frame(family=factor(c(1,2,3,4,5,6,7,9,11,12,13,16,17,18,19,20)),TadSVL=14),interval="confidence",level=0.95,type="percentile",se="boot",bsmethod="wxy",R=1000))
###Data from Table S2, eliminating family = 15 because of only a single ###tadpole.
SirePulseRate <- c(25.74,25.66,26.89,23.94,25.99,26.14,26.11,28.25,25.49,25.66,25.54,26.84,26.78,26.01,25.87,25.69)
TadpoleSurvival <- c(0.68,0.74,0.70,0.26,0.86,0.80,0.62,0.30,0.78,0.82,0.62,0.96,0.58,0.18,0.86,0.86)
Family <- c(1,2,3,4,5,6,7,9,11,12,13,16,17,18,19,20)
FemaleSVL <-  c(44.8,49.7,44.29,48.17,49.72,47.56,43.38,47.96,45.1,48.34,45.72,44.83,47.33,49.38,44.89,47.53)
summary.qr.hof <- data.frame(cbind(Family,SirePulseRate,TadpoleSurvival,FemaleSVL,PredMass.SVL.14.50 = PredMass.SVL.14.qr50[,1],PredMass.SVL.14.50.lwr95 = PredMass.SVL.14.qr50[,2],PredMass.SVL.14.50.upr95 = PredMass.SVL.14.qr50[,3],PredMass.SVL.12.50 = PredMass.SVL.12.qr50[,1],PredMass.SVL.12.50.lwr95 = PredMass.SVL.12.qr50[,2],PredMass.SVL.12.50.upr95 = PredMass.SVL.12.qr50[,3],PredMass.SVL.16.50 = PredMass.SVL.16.qr50[,1],PredMass.SVL.16.50.lwr95 = PredMass.SVL.16.qr50[,2],PredMass.SVL.16.50.upr95 = PredMass.SVL.16.qr50[,3]))
library(Hmisc)
###Median estimates and 95% CI
xYplot(Cbind(PredMass.SVL.14.50,PredMass.SVL.14.50.lwr95,PredMass.SVL.14.50.upr95)~SirePulseRate, data=summary.qr.hof,col=c("black"),ylim=c(0.15,0.35),pch=16,ylab=list(label="Predicted Mass (g) at SVL = 14 mm",cex=1.5),xlab=list(label="Sire Pulse Rate",cex=1.5))
xYplot(Cbind(PredMass.SVL.12.50,PredMass.SVL.12.50.lwr95,PredMass.SVL.12.50.upr95)~SirePulseRate, data=summary.qr.hof,col=c("black"),ylim=c(0.10,0.25),pch=16,ylab="Predicted Mass (g) at SVL = 12 mm",xlab="Sire Pulse Rate",main="Predicted median mass and 95% CI")
xYplot(Cbind(PredMass.SVL.16.50,PredMass.SVL.16.50.lwr95,PredMass.SVL.16.50.upr95)~SirePulseRate, data=summary.qr.hof,col=c("black"),ylim=c(0.20,0.60),pch=16,ylab=list(label="Predicted Mass (g) at SVL = 16 mm",cex=1.5),xlab=list(label="Sire Pulse Rate",cex=1.5),main=list(label="Predicted median mass and 95% CI",cex=1.5))
summary(lm(PredMass.SVL.14.50 ~ SirePulseRate,data=summary.qr.hof))
###Looking at relationship of mass at SVL and survival
xYplot(Cbind(PredMass.SVL.14.50,PredMass.SVL.14.50.lwr95,PredMass.SVL.14.50.upr95)~TadpoleSurvival, data=summary.qr.hof,col=c("black"),ylim=c(0.15,0.35),pch=16,ylab=list(label="Predicted Mass (g) at SVL = 14 mm",cex=1.5),xlab=list(label="Tadpole survival rate",cex=1.5))
summary(lm(PredMass.SVL.14.50 ~ TadpoleSurvival,data=summary.qr.hof))
###Looking at relationship of mass at SVL and Female SVL
xYplot(Cbind(PredMass.SVL.14.50,PredMass.SVL.14.50.lwr95,PredMass.SVL.14.50.upr95)~FemaleSVL, data=summary.qr.hof,col=c("black"),ylim=c(0.15,0.40),pch=16,ylab=list(label="Predicted Mass (g) at SVL = 14 mm",cex=1.5),xlab=list(label="Female SVL (mm)",cex=1.5),main=list(label="Predicted median mass and 95% CI",cex=1.5))
summary(lm(PredMass.SVL.14.50 ~ FemaleSVL,data=summary.qr.hof))