####-----------------------####
#prerequisites
#rm(list=ls())

library(tidyverse) 
library(smatr) #group-specific analysis of regressions, used for species effects in allometric relationship between the lobe structures.
library(lme4) #GLMMs
library(ggResidpanel) #alternative diagnostics tools in std lms
library(car)
library(DHARMa) #GLMM diagnostics.
library(lmerTest) #overrides the original lmer results and supplies p-values
library(ggbeeswarm) #provides dotplots in semi-random distribution.
library(bbmle)
library(boot) # Bootstrapping to estimate confidence intervals for R-squared
library(plotly) #PCA plotting
library(ggfortify)#PCA plotting
library(ggpubr)
library(emmeans) #pairwise species difference as posthoc test.


####-----------------------####
#GLMM and #######                 LMs FIGs 2 and S.
#Model for sex effects

setwd("C:/R-Analysis")

lmdata=read.table(file="Dataset_LOB-FIN-STR_lm-str.txt", sep="\t", header=TRUE)
str(lmdata)

#1. check distribution of volume (all others are factors).
hist(lmdata$volume)
shapiro.test(lmdata$volume)
#it looks relatively normally distributed but shapiro test indicates its not. You can see such a specific species effect and each volume is assigned to
#a specific structure that its not surprising we don't see a completely normal distribution. Z transformation won't change this,
#so we accepted it.

#2. General: Models.
#we have three different sub data sets. Because we don't have age data for all species, and because we don't have sex data for each individual, we tested age and sex effects
#in the appropriate sub data sets, to then make subsequent decisions.

#3. General: Random factors. It is important to note that age, sex, species and clade are in between-subjects predictors, only occuring once in each individual. 
#however structure is an within-subjects predictor, each individual having all 6 categories of the structure predictor in it. so in theory we NEED to account 
#for random slopes. However these models did not converge, likely because we did not have sufficient amounts of single data points. 
#we included the random intercept of individuals, where appropriate.

#4. Sex effects. 
#first test whether sex has an effect as we have quite a limited dataset here. approx 50%.
#first remove all the ones that do not have data for sex.
sex.lmdata=filter(lmdata, sex == "female" | sex == "male")
sex.lmdata

#5. main model to run is following, with structure and sex as test predictors
#, with expected sex-dependent volume differences in different structures.
sex.model=lmer(volume~structure*sex+(1|ID), data=sex.lmdata, REML=F) 
summary(sex.model)

#However, to identify the model with the best fit, we performed model selection on a set of nested models
#, all with the basic structure, all biologically sound.
#lets make a model comparison and find a better model with better diagnostics. 
sex.model2=lmer(volume~structure*sex+spec+(1|ID), data=sex.lmdata, REML=F) 
sex.model3=lmer(volume~structure*sex+age+(1|ID), data=sex.lmdata, REML=F) 
sex.model4=lmer(volume~structure*sex+spec+age+(1|ID), data=sex.lmdata, REML=F) 
sex.model5=lmer(volume~structure*sex+spec*age+(1|ID), data=sex.lmdata, REML=F) 
#most of these models show a very reduced contribution of the random factor, which makes returns a singular/overfitted model.
#AIC comparison (not shown) returns that model1 has a relatively high AIC, though. 
#Hence, assuming that controlling for indvidual differences in the sex differences has a small effect, 
#we continued with linear models/multiple regressions.

sex.model_=lm(volume~structure*sex, data=sex.lmdata) 
sex.model2_=lm(volume~structure*sex+spec, data=sex.lmdata) 
sex.model3_=lm(volume~structure*sex+age, data=sex.lmdata) 
sex.model4_=lm(volume~structure*sex+spec+age, data=sex.lmdata) 
sex.model5_=lm(volume~structure*sex+spec*age, data=sex.lmdata)

AICtab(sex.model_,sex.model2_,sex.model3_,sex.model4_,sex.model5_)

#sex.model5 has the lowest AIC as expected, including the interaction of species and age, but model 4 is minutely higher, but simpler 
#model 4 preferred.

#diagnostics.
resid_panel(sex.model4_, smoother=TRUE)
#indicates largely normal distribution of residuals.
resid_xpanel(sex.model4_)
#the residuals plot equally across the predictors.

summary(sex.model4_)
#this reveals no significant sex differences per structure and overall.

#run the null model.
sex.model_null=lm(volume~structure+spec+age,data=sex.lmdata)
anova(sex.model4_,sex.model_null)
#insignificant.

##sex as factor can be ignored test influential variable.

#extract values.
round(summary(sex.model4_)$fstatistic, digits=3)
round(summary(sex.model4_)$adj.r.squared, digits=3)
coef.sex.effects=round(summary(sex.model4_)$coefficients[ , c(1, 3,4)], digits=3)
coef.sex.effects
write.table(coef.sex.effects,file="Sex.Effects.Coeff.txt",row.names=T,col.names=T,sep="\t")

###

#Model for age effects:
#as other data from the lab (https://doi.org/10.1101/2023.07.11.548546) indicated species specific age effects,
#we tested a three way interaction of structure*age ("structure specific age effects...")*spec ("...between species").

#subset data
age.lmdata=lmdata[lmdata$spec %in% c('erato', 'iulia'), ]
age.lmdata

age.model=lmer(volume~structure*age*spec+(1|ID), data=age.lmdata, REML=F) 
#no sensical alternative nested models available (there is only sex that can be added as a term). 
#Hence, directly to diagnostics.

#alternative nested models.I dont see any after we now know that sex doesnt add much really. 
#diagnostics. 

simulationOutput.age.model <- simulateResiduals(fittedModel = age.model)
plot(simulationOutput.age.model)
#no violation.
#residual distribution alongside (test) predictors.
plotResiduals(simulationOutput.age.model, form = age.lmdata$age)
plotResiduals(simulationOutput.age.model, form = age.lmdata$structure)
plotResiduals(simulationOutput.age.model, form = age.lmdata$spec)
#all good.

#inference
summary(age.model)
round(summary(age.model)$logLik, digits=3)

#lets see whether age-specific species effects (age*spec) are significant in explaining volume variation.
age.null.model.1=lmer(volume~structure*spec+age+(1|ID), data=age.lmdata, REML=F) 
round(summary(age.null.model.1)$logLik, digits=3)
anova(age.null.model.1,age.model,test="Chisq")
#age*spec is a highly significant factor.

#does age itself have an effect?
age.null.model.2=lmer(volume~structure*spec+(1|ID), data=age.lmdata, REML=F) 
round(summary(age.null.model.2)$logLik, digits=3)
anova(age.null.model.2,age.model,test="Chisq")
#yes.

#lets subset this to tease apart these results for each structure as a posthoc test.
age.lmdata.alpha=subset(age.lmdata, structure=="log10_alpha")
age.lmdata.alpha_prime=subset(age.lmdata, structure=="log10_alpha_prime")
age.lmdata.beta=subset(age.lmdata, structure=="log10_beta")
age.lmdata.beta_prime=subset(age.lmdata, structure=="log10_beta_prime")
age.lmdata.gamma=subset(age.lmdata, structure=="log10_gamma")
age.lmdata.Y.lobe=subset(age.lmdata, structure=="log10_Y-lobe")

#models. the variation between individuals is generally now. 
#with these smaller signle structure data sets it makes sense to omit random intercept testing.
#the structure info is still in the variable volume. hence response is "volume".
#models. the variation in ids is so low we can omit it and use standard lms.

age.model.alpha=lm(volume~spec*age,data=age.lmdata.alpha)
summary(age.model.alpha)
coef.age.effects.alpha=summary(age.model.alpha)$coefficients[ , c(1, 3,4)]
coef.age.effects.alpha
#no age effects or spec*age effects.
age.model.alpha_prime=lm(volume~spec*age,data=age.lmdata.alpha_prime)
summary(age.model.alpha_prime)
coef.age.effects.alpha_prime=summary(age.model.alpha_prime)$coefficients[ , c(1, 3,4)]
coef.age.effects.alpha_prime
#age effect.
age.model.beta=lm(volume~spec*age,data=age.lmdata.beta)
summary(age.model.beta)
coef.age.effects.beta=summary(age.model.beta)$coefficients[ , c(1, 3,4)]
coef.age.effects.beta
#no effects.
age.model.beta_prime=lm(volume~spec*age,data=age.lmdata.beta_prime)
summary(age.model.beta_prime)
coef.age.effects.beta_prime=summary(age.model.beta_prime)$coefficients[ , c(1, 3,4)]
coef.age.effects.beta_prime
#age*spec effects!
age.model.gamma=lm(volume~spec*age,data=age.lmdata.gamma)
summary(age.model.gamma)
coef.age.effects.gamma=summary(age.model.gamma)$coefficients[ , c(1, 3,4)]
coef.age.effects.gamma
#age effects.
age.model.Y.lobe=lm(volume~spec*age,data=age.lmdata.Y.lobe)
summary(age.model.Y.lobe)
coef.age.effects.Y.lobe=summary(age.model.Y.lobe)$coefficients[ , c(1, 3,4)]
coef.age.effects.Y.lobe

#age*spec effects.


#plot
age.lmdata$age=factor(age.lmdata$age,levels=
                          c("young","old"))
age.lmdata$spec=factor(age.lmdata$spec,levels=
                           c("iulia","erato"))

age.effects.plot=ggplot(age.lmdata, aes(x=spec, y=volume, fill=age))
age.effects.plot
age.effects.plot2=age.effects.plot+facet_wrap(~structure, scale="free")+
  theme_minimal()
age.effects.plot2
age.effects.plot3=age.effects.plot2+geom_quasirandom(dodge.width=0.75,size=1.5,stroke=0.75,shape=21,aes(colour=spec))
age.effects.plot3
age.effects.plot4=age.effects.plot3+scale_fill_manual(values = c("grey70", "grey20"))
age.effects.plot4
age.effects.plot5=age.effects.plot4+scale_color_manual(values = c("#1293b9ff", "#c22f20ca"))
age.effects.plot5
age.effects.plot6=age.effects.plot5+
  stat_summary(fun.data=mean_sdl, fun.args = list(mult=1), 
                   geom="pointrange", color="black",alpha=0.5, size=1, stroke=0,aes(group=age),position=position_dodge(width=0.75))+
  theme(legend.position = "none")
age.effects.plot6

###

#species effects
#in the young samples only.otherwise we'll bias it against the age effects.

young.lmdata=filter(lmdata, age=="young")
young.lmdata

spec.model=lmer(volume~structure*spec+(1|ID), data=young.lmdata, REML=F) 
round(summary(spec.model)$logLik, digits=3)
#this is our main model. we removed age effects from the data set beforehand, sex doesnt have any effect, 
#and structure*spec is indeed the test variable. so not alternative models for model selection necessary.

#diagnostics. 

simulationOutput.spec.model <- simulateResiduals(fittedModel = spec.model)
plot(simulationOutput.spec.model)
#look good.

#residual distribution alongside (test) predictors.
plotResiduals(simulationOutput.spec.model, form = young.lmdata$structure)
plotResiduals(simulationOutput.spec.model, form = young.lmdata$spec)
#all good.

#can we substitute species with clade, i.e. if there are hardly any differences between Dione juno 
#and Dryas iulia, as well as the two Heliconius species, Clade will offer as much explanation for volume differences
#as species does.

clade.model=lmer(volume~structure*clade+(1|ID), data=young.lmdata, REML=F)
round(summary(clade.model)$logLik, digits=3)
anova(clade.model,spec.model,test="Chisq")

#signficant differences, hence we will want to stay with species.
#lets investigate general effects of species on volume.

spec.null.model=lmer(volume~structure+spec+(1|ID), data=young.lmdata, REML=F) 
round(summary(spec.null.model)$logLik, digits=3)
anova(spec.null.model,spec.model,test="Chisq")
round(anova(spec.null.model,spec.model,test="Chisq")$Chisq, digits=3)
#very highly significant effects of the interaction of structure and species.


#posthoc structure specific tests

posthocspecies.structure <- emmeans(spec.model, list(pairwise ~ spec|structure),adjust="tukey")
plot(posthocspecies.structure)
emmeans=as.data.frame(summary(posthocspecies.structure$`pairwise differences of spec | structure`))


#plot
young.lmdata$spec=factor(young.lmdata$spec,levels=
                           c("iulia","juno","erato","melpomene"))

spec.effects.plotA=ggplot(young.lmdata, aes(x=structure, y=volume, fill=spec)) + 
  theme_minimal()
spec.effects.plotA
spec.effects.plot2A=spec.effects.plotA+geom_beeswarm(dodge.width=0.88,size=1.3,aes(colour=spec))
spec.effects.plot2A
spec.effects.plot3A=spec.effects.plot2A+scale_colour_manual(values = c("#1293b9ff", "#7dadb5ff","#c22f20ca","#cb8781ff"))+
  theme(legend.position="none")
spec.effects.plot3A
spec.effects.plot4A=spec.effects.plot3A+
  stat_summary(fun.data=mean_sdl, fun.args = list(mult=1), 
               geom="pointrange", color="black",alpha=0.5, size=1, stroke=0,aes(group=spec),position=position_dodge(width=0.88))+
  theme(legend.position = "none")
spec.effects.plot4A

#effect sizes of single structure lms. R2 is retreivable from the lms of course, but how do we get confidence intervals
#of them? we used a bootstrapping approach, supported by ChatGPT.

#first subset data and run models
young.lmdata.alpha=subset(young.lmdata, structure=="log10_alpha")
young.lmdata.alpha_prime=subset(young.lmdata, structure=="log10_alpha_prime")
young.lmdata.beta=subset(young.lmdata, structure=="log10_beta")
young.lmdata.beta_prime=subset(young.lmdata, structure=="log10_beta_prime")
young.lmdata.gamma=subset(young.lmdata, structure=="log10_gamma")
young.lmdata.Y.lobe=subset(young.lmdata, structure=="log10_Y-lobe")

spec.model.alpha=lm(volume~spec,data=young.lmdata.alpha)
summary(spec.model.alpha)
spec.model.alpha_prime=lm(volume~spec,data=young.lmdata.alpha_prime)
summary(spec.model.alpha_prime)
spec.model.beta=lm(volume~spec,data=young.lmdata.beta)
summary(spec.model.beta)
spec.model.beta_prime=lm(volume~spec,data=young.lmdata.beta_prime)
summary(spec.model.beta_prime)
spec.model.gamma=lm(volume~spec,data=young.lmdata.gamma)
summary(spec.model.gamma)
spec.model.Y.lobe=lm(volume~spec,data=young.lmdata.Y.lobe)
summary(spec.model.Y.lobe)

# Function to calculate R-squared
calc_r_squared <- function(data, indices) {
  fit <- lm(volume ~ spec, data = data[indices,])
  return(summary(fit)$r.squared)
}

# Bootstrapping to estimate confidence intervals for R-squared
boot_results.alpha <- boot(data = young.lmdata.alpha, statistic = calc_r_squared, R = 1000)

# 95% confidence intervals for R-squared
conf_int_r_squared.alpha <- boot.ci(boot_results.alpha, type = "bca")
print(conf_int_r_squared.alpha)

#for alpha-prime
boot_results.alpha_prime <- boot(data = young.lmdata.m.alpha_prime, statistic = calc_r_squared, R = 1000)

# 95% confidence intervals for R-squared
conf_int_r_squared.alpha_prime <- boot.ci(boot_results.alpha_prime, type = "bca")
print(conf_int_r_squared.alpha_prime)

#for beta
boot_results.beta <- boot(data = young.lmdata.m.beta, statistic = calc_r_squared, R = 1000)

# 95% confidence intervals for R-squared
conf_int_r_squared.beta <- boot.ci(boot_results.beta, type = "bca")
print(conf_int_r_squared.beta)

#for beta prime
boot_results.beta_prime <- boot(data = young.lmdata.m.beta_prime, statistic = calc_r_squared, R = 1000)

# 95% confidence intervals for R-squared
conf_int_r_squared.beta_prime <- boot.ci(boot_results.beta_prime, type = "bca")
print(conf_int_r_squared.beta_prime)

#for gamma
boot_results.gamma <- boot(data = young.lmdata.m.gamma, statistic = calc_r_squared, R = 1000)

# 95% confidence intervals for R-squared
conf_int_r_squared.gamma <- boot.ci(boot_results.gamma, type = "bca")
print(conf_int_r_squared.gamma)

#for Y lobe
boot_results.Y.lobe <- boot(data = young.lmdata.m.Y.lobe, statistic = calc_r_squared, R = 1000)

# 95% confidence intervals for R-squared
conf_int_r_squared.Y.lobe <- boot.ci(boot_results.Y.lobe, type = "bca")
print(conf_int_r_squared.Y.lobe)

#manual export of data.

Rsq_CIs_data=read.table(file="Rsq_CIs.txt", sep="\t", header=TRUE)
str(Rsq_CIs_data)

Rsq_CIs_data$structure=factor(Rsq_CIs_data$structure,levels=
                                c("alpha","beta_prime","gamma","alpha_prime","beta","Y.lobe"))
Rsq_CI.plot= ggplot(Rsq_CIs_data, aes(structure, rsquared)) +
  geom_point(aes(colour = structure),size=4.5) + theme_minimal()+
  geom_errorbar(aes(ymin = lowerCI, ymax = higherCI,colour = structure),width=0.3)
Rsq_CI.plot
Rsq_CI.plot2=Rsq_CI.plot+scale_colour_manual(values=c("#5fd2f2ff","#ffe646ff","#2e98c0ff","#ddc657ff","#ff6d0cff","#8033aaff"))+
  theme(legend.position = "none")
Rsq_CI.plot2

####-----------------------####
#PCAs
#library(plotly)
#library(ggfortify)

sdata=read.table(file="Dataset_LOB-FIN-STR_.txt", sep="\t", header=TRUE)
str(sdata)

#filter for only young ones (removing species bias)
sdata.young=filter(sdata, age == "young")
str(sdata.young)

#now remove the non-numerical categories.
lob.data=sdata.young[, 5:10]
str(lob.data)

#PCA
pca1=prcomp(lob.data)
summary(pca1)#prints summary
pca1.importance=summary(pca1)$importance
pca1.values=pca1$rotation

##plot result. pc1 and 2.
pca1.plot=autoplot(pca1,data=sdata.young, colour="spec",size=4.5, alpha=0.7, 
                   loadings = TRUE, loadings.colour = 'darkgrey',
                   loadings.label = TRUE, loadings.label.size = 4)+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),legend.position = "none")+
  scale_colour_manual(values = c("#c22f20ca", "#1293b9ff", "#7dadb5ff","#cb8781ff"))
pca1.plot

####-----------------------####
#SMATr ANALYSIS
#selection for young ones only. so use sdata.young
str(sdata.young)

#SMATr test itself.

Elev.test1=sma(log10_gamma~log10_alpha_prime+spec,multcomp=T, robust=T, data=sdata.young)
summary(Elev.test1)
Elev.test1.grpsum=round(Elev.test1$groupsummary[4], digits=3)
Elev.test1.multcomp=Elev.test1$multcompresult[ , c(3,4)]
Elev.test1$gtr #to extract Wald stat.

Elev.test2=sma(log10_gamma~log10_beta_prime+spec,multcomp=T,robust=T,data=sdata.young)
summary(Elev.test2)
Elev.test2.grpsum=round(Elev.test2$groupsummary[4], digits=3)
Elev.test2.multcomp=Elev.test2$multcompresult[ , c(3,4)]
Elev.test2$gtr

Elev.test3=sma(log10_gamma~log10_alpha+spec,multcomp=T,robust=T,data=sdata.young)
summary(Elev.test3)
Elev.test3.grpsum=round(Elev.test3$groupsummary[4], digits=3)
Elev.test3.multcomp=Elev.test3$multcompresult[ , c(3,4)]
Elev.test3$gtr

Elev.test4=sma(log10_gamma~log10_beta+spec,multcomp=T,robust=T,data=sdata.young)
summary(Elev.test4)
Elev.test4.grpsum=round(Elev.test4$groupsummary[4], digits=3)
Elev.test4.multcomp=Elev.test4$multcompresult[ , c(3,4)]
Elev.test4$gtr

Elev.test5=sma(log10_gamma~log10_Y.lobe+spec,multcomp=T,robust=T,data=sdata.young)
summary(Elev.test5)
Elev.test5.grpsum=round(Elev.test5$groupsummary[4], digits=3)
Elev.test5.multcomp=Elev.test5$multcompresult[ , c(3,4)]
Elev.test5$gtr

Elev.test6=sma(log10_alpha_prime~log10_gamma+spec,multcomp=T,robust=T,data=sdata.young)
summary(Elev.test6)
Elev.test6.grpsum=round(Elev.test6$groupsummary[4], digits=3)
Elev.test6.multcomp=Elev.test6$multcompresult[ , c(3,4)]
Elev.test6$gtr

Elev.test7=sma(log10_alpha_prime~log10_beta_prime+spec,multcomp=T,robust=T,data=sdata.young)
summary(Elev.test7)
Elev.test7.grpsum=round(Elev.test7$groupsummary[4], digits=3)
Elev.test7.multcomp=Elev.test7$multcompresult[ , c(3,4)]
Elev.test7$gtr

Elev.test8=sma(log10_alpha_prime~log10_alpha+spec,multcomp=T,robust=T,data=sdata.young)
summary(Elev.test8)
Elev.test8.grpsum=round(Elev.test8$groupsummary[4], digits=3)
Elev.test8.multcomp=Elev.test8$multcompresult[ , c(3,4)]
Elev.test8$gtr

Elev.test9=sma(log10_alpha_prime~log10_beta+spec,multcomp=T,robust=T,data=sdata.young)
summary(Elev.test9)
Elev.test9.grpsum=round(Elev.test9$groupsummary[4], digits=3)
Elev.test9.multcomp=Elev.test9$multcompresult[ , c(3,4)]
Elev.test9$gtr

Elev.test10=sma(log10_alpha_prime~log10_Y.lobe+spec,multcomp=T,robust=T,data=sdata.young)
summary(Elev.test10)
Elev.test10.grpsum=round(Elev.test10$groupsummary[4], digits=3)
Elev.test10.multcomp=Elev.test10$multcompresult[ , c(3,4)]
Elev.test10$gtr

Elev.test11=sma(log10_beta_prime~log10_gamma+spec,multcomp=T,robust=T,data=sdata.young)
summary(Elev.test11)
Elev.test11.grpsum=round(Elev.test11$groupsummary[4], digits=3)
Elev.test11.multcomp=Elev.test11$multcompresult[ , c(3,4)]
Elev.test11$gtr

Elev.test12=sma(log10_beta_prime~log10_alpha_prime+spec,multcomp=T,robust=T,data=sdata.young)
summary(Elev.test12)
Elev.test12.grpsum=round(Elev.test12$groupsummary[4], digits=3)
Elev.test12.multcomp=Elev.test12$multcompresult[ , c(3,4)]
Elev.test12$gtr

Elev.test13=sma(log10_beta_prime~log10_alpha+spec,multcomp=T,robust=T,data=sdata.young)
summary(Elev.test13)
Elev.test13.grpsum=round(Elev.test13$groupsummary[4], digits=3)
Elev.test13.multcomp=Elev.test13$multcompresult[ , c(3,4)]
Elev.test13$gtr

Elev.test14=sma(log10_beta_prime~log10_beta+spec,multcomp=T,robust=T,data=sdata.young)
summary(Elev.test14)
Elev.test14.grpsum=round(Elev.test14$groupsummary[4], digits=3)
Elev.test14.multcomp=Elev.test14$multcompresult[ , c(3,4)]
Elev.test14$gtr

Elev.test15=sma(log10_beta_prime~log10_Y.lobe+spec,multcomp=T,robust=T,data=sdata.young)
summary(Elev.test15)
Elev.test15.grpsum=round(Elev.test15$groupsummary[4], digits=3)
Elev.test15.multcomp=Elev.test15$multcompresult[ , c(3,4)]
Elev.test15$gtr

Elev.test16=sma(log10_alpha~log10_gamma+spec,multcomp=T,robust=T,data=sdata.young)
summary(Elev.test16)
Elev.test16.grpsum=round(Elev.test16$groupsummary[4], digits=3)
Elev.test16.multcomp=Elev.test16$multcompresult[ , c(3,4)]
Elev.test16$gtr

Elev.test17=sma(log10_alpha~log10_alpha_prime+spec,multcomp=T,robust=T,data=sdata.young)
summary(Elev.test17)
Elev.test17.grpsum=round(Elev.test17$groupsummary[4], digits=3)
Elev.test17.multcomp=Elev.test17$multcompresult[ , c(3,4)]
Elev.test17$gtr

Elev.test18=sma(log10_alpha~log10_beta_prime+spec,multcomp=T,robust=T,data=sdata.young)
summary(Elev.test18)
Elev.test18.grpsum=round(Elev.test18$groupsummary[4], digits=3)
Elev.test18.multcomp=Elev.test18$multcompresult[ , c(3,4)]
Elev.test18$gtr

Elev.test19=sma(log10_alpha~log10_beta+spec,multcomp=T,robust=T,data=sdata.young)
summary(Elev.test19)
Elev.test19.grpsum=round(Elev.test19$groupsummary[4], digits=3)
Elev.test19.multcomp=Elev.test19$multcompresult[ , c(3,4)]
Elev.test19$gtr

Elev.test20=sma(log10_alpha~log10_Y.lobe+spec,multcomp=T,robust=T,data=sdata.young)
summary(Elev.test20)
Elev.test20.grpsum=round(Elev.test20$groupsummary[4], digits=3)
Elev.test20.multcomp=Elev.test20$multcompresult[ , c(3,4)]
Elev.test20$gtr

Elev.test21=sma(log10_beta~log10_gamma+spec,multcomp=T,robust=T,data=sdata.young)
summary(Elev.test21)
Elev.test21.grpsum=round(Elev.test21$groupsummary[4], digits=3)
Elev.test21.multcomp=Elev.test21$multcompresult[ , c(3,4)]
Elev.test21$gtr

Elev.test22=sma(log10_beta~log10_alpha_prime+spec,multcomp=T,robust=T,data=sdata.young)
summary(Elev.test22)
Elev.test22.grpsum=round(Elev.test22$groupsummary[4], digits=3)
Elev.test22.multcomp=Elev.test22$multcompresult[ , c(3,4)]
Elev.test22$gtr

Elev.test23=sma(log10_beta~log10_beta_prime+spec,multcomp=T,robust=T,data=sdata.young)
summary(Elev.test23)
Elev.test23.grpsum=round(Elev.test23$groupsummary[4], digits=3)
Elev.test23.multcomp=Elev.test23$multcompresult[ , c(3,4)]
Elev.test23$gtr

Elev.test24=sma(log10_beta~log10_alpha+spec,multcomp=T,robust=T,data=sdata.young)
summary(Elev.test24)
Elev.test24.grpsum=round(Elev.test24$groupsummary[4], digits=3)
Elev.test24.multcomp=Elev.test24$multcompresult[ , c(3,4)]
Elev.test24$gtr

Elev.test25=sma(log10_beta~log10_Y.lobe+spec,multcomp=T,robust=T,data=sdata.young)
summary(Elev.test25)
Elev.test25.grpsum=round(Elev.test25$groupsummary[4], digits=3)
Elev.test25.multcomp=Elev.test25$multcompresult[ , c(3,4)]
Elev.test25$gtr

Elev.test26=sma(log10_Y.lobe~log10_gamma+spec,multcomp=T,robust=T,data=sdata.young)
summary(Elev.test26)
Elev.test26.grpsum=round(Elev.test26$groupsummary[4], digits=3)
Elev.test26.multcomp=Elev.test26$multcompresult[ , c(3,4)]
Elev.test26$gtr

Elev.test27=sma(log10_Y.lobe~log10_alpha_prime+spec,multcomp=T,robust=T,data=sdata.young)
summary(Elev.test27)
Elev.test27.grpsum=round(Elev.test27$groupsummary[4], digits=3)
Elev.test27.multcomp=Elev.test27$multcompresult[ , c(3,4)]
Elev.test27$gtr

Elev.test28=sma(log10_Y.lobe~log10_beta_prime+spec,multcomp=T,robust=T,data=sdata.young)
summary(Elev.test28)
Elev.test28.grpsum=round(Elev.test28$groupsummary[4], digits=3)
Elev.test28.multcomp=Elev.test28$multcompresult[ , c(3,4)]
Elev.test28$gtr

Elev.test29=sma(log10_Y.lobe~log10_alpha+spec,multcomp=T,robust=T,data=sdata.young)
summary(Elev.test29)
Elev.test29.grpsum=round(Elev.test29$groupsummary[4], digits=3)
Elev.test29.multcomp=Elev.test29$multcompresult[ , c(3,4)]
Elev.test29$gtr

Elev.test30=sma(log10_Y.lobe~log10_beta+spec,multcomp=T,robust=T,data=sdata.young)
summary(Elev.test30)
Elev.test30.grpsum=round(Elev.test30$groupsummary[4], digits=3)
Elev.test30.multcomp=Elev.test30$multcompresult[ , c(3,4)]
Elev.test30$gtr

#add the lines to the plots.

sma1.coefs=data.frame(coef(Elev.test1))
sma1.coefs
sma2.coefs=data.frame(coef(Elev.test2))
sma2.coefs
sma3.coefs=data.frame(coef(Elev.test3))
sma3.coefs
sma4.coefs=data.frame(coef(Elev.test4))
sma4.coefs
sma5.coefs=data.frame(coef(Elev.test5))
sma5.coefs

sma6.coefs=data.frame(coef(Elev.test6))
sma6.coefs
sma7.coefs=data.frame(coef(Elev.test7))
sma7.coefs
sma8.coefs=data.frame(coef(Elev.test8))
sma8.coefs
sma9.coefs=data.frame(coef(Elev.test9))
sma9.coefs
sma10.coefs=data.frame(coef(Elev.test10))
sma10.coefs

sma11.coefs=data.frame(coef(Elev.test11))
sma11.coefs
sma12.coefs=data.frame(coef(Elev.test12))
sma12.coefs
sma13.coefs=data.frame(coef(Elev.test13))
sma13.coefs
sma14.coefs=data.frame(coef(Elev.test14))
sma14.coefs
sma15.coefs=data.frame(coef(Elev.test15))
sma15.coefs

sma16.coefs=data.frame(coef(Elev.test16))
sma16.coefs
sma17.coefs=data.frame(coef(Elev.test17))
sma17.coefs
sma18.coefs=data.frame(coef(Elev.test18))
sma18.coefs
sma19.coefs=data.frame(coef(Elev.test19))
sma19.coefs
sma20.coefs=data.frame(coef(Elev.test20))
sma20.coefs

sma21.coefs=data.frame(coef(Elev.test21))
sma21.coefs
sma22.coefs=data.frame(coef(Elev.test22))
sma22.coefs
sma23.coefs=data.frame(coef(Elev.test23))
sma23.coefs
sma24.coefs=data.frame(coef(Elev.test24))
sma24.coefs
sma25.coefs=data.frame(coef(Elev.test25))
sma25.coefs

sma26.coefs=data.frame(coef(Elev.test26))
sma26.coefs
sma27.coefs=data.frame(coef(Elev.test27))
sma27.coefs
sma28.coefs=data.frame(coef(Elev.test28))
sma28.coefs
sma29.coefs=data.frame(coef(Elev.test29))
sma29.coefs
sma30.coefs=data.frame(coef(Elev.test30))
sma30.coefs

lobp1.plot=ggplot(sdata.young,aes(x=log10_alpha_prime, y=log10_gamma))+geom_point(aes(col=spec,size=0.04, alpha=0.5))+
  xlab(expression(bold(log[10](alpha~"'"))))+ylab(expression(bold(log[10](gamma))))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.border = element_blank(), panel.background = element_blank(),axis.line = element_line(color = 'black',size=0.8), legend.position = "none")
lobp1.plot
lobp2.plot=ggplot(sdata.young,aes(x=log10_beta_prime, y=log10_gamma))+geom_point(aes(col=spec, size=0.04, alpha=0.5))+
  xlab(expression(bold(log[10](beta~"'"))))+ylab(expression(bold(log[10](gamma))))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.border = element_blank(), panel.background = element_blank(),axis.line = element_line(color = 'black',size=0.8), legend.position = "none")
lobp2.plot
lobp3.plot=ggplot(sdata.young,aes(x=log10_alpha, y=log10_gamma))+geom_point(aes(col=spec,  size=0.04, alpha=0.5))+
  xlab(expression(bold(log[10](alpha))))+ylab(expression(bold(log[10](gamma))))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.border = element_blank(), panel.background = element_blank(),axis.line = element_line(color = 'black',size=0.8), legend.position = "none")
lobp3.plot
lobp4.plot=ggplot(sdata.young,aes(x=log10_beta, y=log10_gamma))+geom_point(aes(col=spec,  size=0.04, alpha=0.5))+
  xlab(expression(bold(log[10](beta))))+ylab(expression(bold(log[10](gamma))))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.border = element_blank(), panel.background = element_blank(),axis.line = element_line(color = 'black',size=0.8), legend.position = "none")
lobp4.plot
lobp5.plot=ggplot(sdata.young,aes(x=log10_Y.lobe, y=log10_gamma))+geom_point(aes(col=spec,size=0.04, alpha=0.5))+
  xlab(expression(bold(log[10](Y))))+ylab(expression(bold(log[10](gamma))))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.border = element_blank(), panel.background = element_blank(),axis.line = element_line(color = 'black',size=0.8), legend.position = "none")
lobp5.plot

lobp6.plot=ggplot(sdata.young,aes(x=log10_gamma, y=log10_alpha_prime))+geom_point(aes(col=spec, size=0.04, alpha=0.5))+
  xlab(expression(bold(log[10](gamma))))+ylab(expression(bold(log[10](alpha~"'"))))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.border = element_blank(), panel.background = element_blank(),axis.line = element_line(color = 'black',size=0.8), legend.position = "none")
lobp6.plot
lobp7.plot=ggplot(sdata.young,aes(x=log10_beta_prime, y=log10_alpha_prime))+geom_point(aes(col=spec,  size=0.04, alpha=0.5))+
  xlab(expression(bold(log[10](beta~"'"))))+ylab(expression(bold(log[10](alpha~"'"))))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.border = element_blank(), panel.background = element_blank(),axis.line = element_line(color = 'black',size=0.8), legend.position = "none")
lobp7.plot
lobp8.plot=ggplot(sdata.young,aes(x=log10_alpha, y=log10_alpha_prime))+geom_point(aes(col=spec, size=0.04, alpha=0.5))+
  xlab(expression(bold(log[10](alpha))))+ylab(expression(bold(log[10](alpha~"'"))))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.border = element_blank(), panel.background = element_blank(),axis.line = element_line(color = 'black',size=0.8), legend.position = "none")
lobp8.plot
lobp9.plot=ggplot(sdata.young,aes(x=log10_beta, y=log10_alpha_prime))+geom_point(aes(col=spec, size=0.04, alpha=0.5))+
  xlab(expression(bold(log[10](beta))))+ylab(expression(bold(log[10](alpha~"'"))))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.border = element_blank(), panel.background = element_blank(),axis.line = element_line(color = 'black',size=0.8), legend.position = "none")
lobp9.plot
lobp10.plot=ggplot(sdata.young,aes(x=log10_Y.lobe, y=log10_alpha_prime))+geom_point(aes(col=spec, size=0.04, alpha=0.5))+
  xlab(expression(bold(log[10](Y))))+ylab(expression(bold(log[10](alpha~"'"))))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.border = element_blank(), panel.background = element_blank(),axis.line = element_line(color = 'black',size=0.8), legend.position = "none")
lobp10.plot


lobp11.plot=ggplot(sdata.young,aes(x=log10_gamma, y=log10_beta_prime))+geom_point(aes(col=spec,size=0.04, alpha=0.5))+
  xlab(expression(bold(log[10](gamma))))+ylab(expression(bold(log[10](beta~"'"))))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.border = element_blank(), panel.background = element_blank(),axis.line = element_line(color = 'black',size=0.8), legend.position = "none")
lobp11.plot
lobp12.plot=ggplot(sdata.young,aes(x=log10_alpha_prime, y=log10_beta_prime))+geom_point(aes(col=spec, size=0.04, alpha=0.5))+
  xlab(expression(bold(log[10](alpha~"'"))))+ylab(expression(bold(log[10](beta~"'"))))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.border = element_blank(), panel.background = element_blank(),axis.line = element_line(color = 'black',size=0.8), legend.position = "none")
lobp12.plot
lobp13.plot=ggplot(sdata.young,aes(x=log10_alpha, y=log10_beta_prime))+geom_point(aes(col=spec, size=0.04, alpha=0.5))+
  xlab(expression(bold(log[10](alpha))))+ylab(expression(bold(log[10](beta~"'"))))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.border = element_blank(), panel.background = element_blank(),axis.line = element_line(color = 'black',size=0.8), legend.position = "none")
lobp13.plot
lobp14.plot=ggplot(sdata.young,aes(x=log10_beta, y=log10_beta_prime))+geom_point(aes(col=spec, size=0.04, alpha=0.5))+
  xlab(expression(bold(log[10](beta))))+ylab(expression(bold(log[10](beta~"'"))))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.border = element_blank(), panel.background = element_blank(),axis.line = element_line(color = 'black',size=0.8), legend.position = "none")
lobp14.plot
lobp15.plot=ggplot(sdata.young,aes(x=log10_Y.lobe, y=log10_beta_prime))+geom_point(aes(col=spec,  size=0.04, alpha=0.5))+
  xlab(expression(bold(log[10](Y))))+ylab(expression(bold(log[10](beta~"'"))))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.border = element_blank(), panel.background = element_blank(),axis.line = element_line(color = 'black',size=0.8), legend.position = "none")
lobp15.plot

lobp16.plot=ggplot(sdata.young,aes(x=log10_gamma, y=log10_alpha))+geom_point(aes(col=spec, size=0.04, alpha=0.5))+
  xlab(expression(bold(log[10](gamma))))+ylab(expression(bold(log[10](alpha))))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.border = element_blank(), panel.background = element_blank(),axis.line = element_line(color = 'black',size=0.8), legend.position = "none")
lobp16.plot
lobp17.plot=ggplot(sdata.young,aes(x=log10_alpha_prime, y=log10_alpha))+geom_point(aes(col=spec,size=0.04, alpha=0.5))+
  xlab(expression(bold(log[10](alpha~"'"))))+ylab(expression(bold(log[10](alpha))))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.border = element_blank(), panel.background = element_blank(),axis.line = element_line(color = 'black',size=0.8), legend.position = "none")
lobp17.plot
lobp18.plot=ggplot(sdata.young,aes(x=log10_beta_prime, y=log10_alpha))+geom_point(aes(col=spec,size=0.04, alpha=0.5))+
  xlab(expression(bold(log[10](beta~"'"))))+ylab(expression(bold(log[10](alpha))))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.border = element_blank(), panel.background = element_blank(),axis.line = element_line(color = 'black',size=0.8), legend.position = "none")
lobp18.plot
lobp19.plot=ggplot(sdata.young,aes(x=log10_beta, y=log10_alpha))+geom_point(aes(col=spec, size=0.04, alpha=0.5))+
  xlab(expression(bold(log[10](beta))))+ylab(expression(bold(log[10](alpha))))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.border = element_blank(), panel.background = element_blank(),axis.line = element_line(color = 'black',size=0.8), legend.position = "none")
lobp19.plot
lobp20.plot=ggplot(sdata.young,aes(x=log10_Y.lobe, y=log10_alpha))+geom_point(aes(col=spec, size=0.04, alpha=0.5))+
  xlab(expression(bold(log[10](Y))))+ylab(expression(bold(log[10](alpha))))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.border = element_blank(), panel.background = element_blank(),axis.line = element_line(color = 'black',size=0.8), legend.position = "none")
lobp20.plot

lobp21.plot=ggplot(sdata.young,aes(x=log10_gamma, y=log10_beta))+geom_point(aes(col=spec,  size=0.04, alpha=0.5))+
  xlab(expression(bold(log[10](gamma))))+ylab(expression(bold(log[10](beta))))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.border = element_blank(), panel.background = element_blank(),axis.line = element_line(color = 'black',size=0.8), legend.position = "none")
lobp21.plot
lobp22.plot=ggplot(sdata.young,aes(x=log10_alpha_prime, y=log10_beta))+geom_point(aes(col=spec, size=0.04, alpha=0.5))+
  xlab(expression(bold(log[10](alpha~"'"))))+ylab(expression(bold(log[10](beta))))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.border = element_blank(), panel.background = element_blank(),axis.line = element_line(color = 'black',size=0.8), legend.position = "none")
lobp22.plot
lobp23.plot=ggplot(sdata.young,aes(x=log10_beta_prime, y=log10_beta))+geom_point(aes(col=spec, size=0.04, alpha=0.5))+
  xlab(expression(bold(log[10](beta~"'"))))+ylab(expression(bold(log[10](beta))))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.border = element_blank(), panel.background = element_blank(),axis.line = element_line(color = 'black',size=0.8), legend.position = "none")
lobp23.plot
lobp24.plot=ggplot(sdata.young,aes(x=log10_alpha, y=log10_beta))+geom_point(aes(col=spec, size=0.04, alpha=0.5))+
  xlab(expression(bold(log[10](alpha))))+ylab(expression(bold(log[10](beta))))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.border = element_blank(), panel.background = element_blank(),axis.line = element_line(color = 'black',size=0.8), legend.position = "none")
lobp24.plot
lobp25.plot=ggplot(sdata.young,aes(x=log10_Y.lobe, y=log10_beta))+geom_point(aes(col=spec, size=0.04, alpha=0.5))+
  xlab(expression(bold(log[10](Y))))+ylab(expression(bold(log[10](beta))))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.border = element_blank(), panel.background = element_blank(),axis.line = element_line(color = 'black',size=0.8), legend.position = "none")
lobp25.plot


lobp26.plot=ggplot(sdata.young,aes(x=log10_gamma, y=log10_Y.lobe))+geom_point(aes(col=spec,size=0.04, alpha=0.5))+
  xlab(expression(bold(log[10](gamma))))+ylab(expression(bold(log[10](Y))))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.border = element_blank(), panel.background = element_blank(),axis.line = element_line(color = 'black',size=0.8), legend.position = "none")
lobp26.plot
lobp27.plot=ggplot(sdata.young,aes(x=log10_alpha_prime, y=log10_Y.lobe))+geom_point(aes(col=spec, size=0.04, alpha=0.5))+
  xlab(expression(bold(log[10](alpha~"'"))))+ylab(expression(bold(log[10](Y))))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.border = element_blank(), panel.background = element_blank(),axis.line = element_line(color = 'black',size=0.8), legend.position = "none")
lobp27.plot
lobp28.plot=ggplot(sdata.young,aes(x=log10_beta_prime, y=log10_Y.lobe))+geom_point(aes(col=spec, size=0.04, alpha=0.5))+
  xlab(expression(bold(log[10](beta~"'"))))+ylab(expression(bold(log[10](Y))))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.border = element_blank(), panel.background = element_blank(),axis.line = element_line(color = 'black',size=0.8), legend.position = "none")
lobp28.plot
lobp29.plot=ggplot(sdata.young,aes(x=log10_alpha, y=log10_Y.lobe))+geom_point(aes(col=spec, size=0.04, alpha=0.5))+
  xlab(expression(bold(log[10](alpha))))+ylab(expression(bold(log[10](Y))))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.border = element_blank(), panel.background = element_blank(),axis.line = element_line(color = 'black',size=0.8), legend.position = "none")
lobp29.plot
lobp30.plot=ggplot(sdata.young,aes(x=log10_beta, y=log10_Y.lobe))+geom_point(aes(col=spec, size=0.04, alpha=0.5))+
  xlab(expression(bold(log[10](beta))))+ylab(expression(bold(log[10](Y))))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.border = element_blank(), panel.background = element_blank(),axis.line = element_line(color = 'black',size=0.8), legend.position = "none")
lobp30.plot

lobp1.1.plot=ggplot(sdata.young,aes(x=log10_gamma, y=log10_gamma))+geom_point(aes(size=0.04, alpha=0.5),color='grey')+
  xlab(expression(bold(log[10](gamma))))+ylab(expression(bold(log[10](gamma))))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.border = element_blank(), panel.background = element_blank(),axis.line = element_line(color = 'black',size=0.8), legend.position = "none")
lobp1.1.plot

lobp8.1.plot=ggplot(sdata.young,aes(x=log10_alpha_prime, y=log10_alpha_prime))+geom_point(aes(size=0.04, alpha=0.5),color='grey')+
  xlab(expression(bold(log[10](alpha~"'"))))+ylab(expression(bold(log[10](alpha~"'"))))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.border = element_blank(), panel.background = element_blank(),axis.line = element_line(color = 'black',size=0.8), legend.position = "none")
lobp8.1.plot

lobp15.1.plot=ggplot(sdata.young,aes(x=log10_beta_prime, y=log10_beta_prime))+geom_point(aes(size=0.04, alpha=0.5),color='grey')+
  xlab(expression(bold(log[10](beta~"'"))))+ylab(expression(bold(log[10](beta~"'"))))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.border = element_blank(), panel.background = element_blank(),axis.line = element_line(color = 'black',size=0.8), legend.position = "none")
lobp15.1.plot

lobp22.1.plot=ggplot(sdata.young,aes(x=log10_alpha, y=log10_alpha))+geom_point(aes(size=0.04, alpha=0.5),color='grey')+
  xlab(expression(bold(log[10](alpha))))+ylab(expression(bold(log[10](alpha))))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.border = element_blank(), panel.background = element_blank(),axis.line = element_line(color = 'black',size=0.8), legend.position = "none")
lobp22.1.plot

lobp29.1.plot=ggplot(sdata.young,aes(x=log10_beta, y=log10_beta))+geom_point(aes(size=0.04, alpha=0.5),color='grey')+
  xlab(expression(bold(log[10](beta))))+ylab(expression(bold(log[10](beta))))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.border = element_blank(), panel.background = element_blank(),axis.line = element_line(color = 'black',size=0.8), legend.position = "none")
lobp29.1.plot

lobp36.1.plot=ggplot(sdata.young,aes(x=log10_Y.lobe, y=log10_Y.lobe))+geom_point(aes( size=0.04, alpha=0.5),color='grey')+
  xlab(expression(bold(log[10](Y))))+ylab(expression(bold(log[10](Y))))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.border = element_blank(), panel.background = element_blank(),axis.line = element_line(color = 'black',size=0.8), legend.position = "none")
lobp36.1.plot


lobp1.plot.l=lobp1.plot+
  geom_segment(aes(x = min(log10_alpha_prime[which(spec=='erato')]), xend = max(log10_alpha_prime[which(spec=='erato')]), 
                   y = sma1.coefs[1,1] + sma1.coefs[1,2]*min(log10_alpha_prime[which(spec=='erato')]), yend = sma1.coefs[1,1]+
                     sma1.coefs[1,2]*max(log10_alpha_prime[which(spec=='erato')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_alpha_prime[which(spec=='iulia')]), xend = max(log10_alpha_prime[which(spec=='iulia')]), 
                   y = sma1.coefs[2,1] + sma1.coefs[2,2]*min(log10_alpha_prime[which(spec=='iulia')]), yend = sma1.coefs[2,1]+
                     sma1.coefs[2,2]*max(log10_alpha_prime[which(spec=='iulia')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_alpha_prime[which(spec=='juno')]), xend = max(log10_alpha_prime[which(spec=='juno')]), 
                   y = sma1.coefs[3,1] + sma1.coefs[3,2]*min(log10_alpha_prime[which(spec=='juno')]), yend = sma1.coefs[3,1]+
                     sma1.coefs[3,2]*max(log10_alpha_prime[which(spec=='juno')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_alpha_prime[which(spec=='melpomene')]), xend = max(log10_alpha_prime[which(spec=='melpomene')]), 
                   y = sma1.coefs[4,1] + sma1.coefs[4,2]*min(log10_alpha_prime[which(spec=='melpomene')]), yend = sma1.coefs[4,1]+
                     sma1.coefs[4,2]*max(log10_alpha_prime[which(spec=='melpomene')])),colour="darkgrey",linewidth=1)+
  scale_colour_manual(values = c("#c22f20ca", "#1293b9ff", "#7dadb5ff","#cb8781ff"))
lobp1.plot.l

lobp2.plot.l=lobp2.plot+
  geom_segment(aes(x = min(log10_beta_prime[which(spec=='erato')]), xend = max(log10_beta_prime[which(spec=='erato')]), 
                   y = sma2.coefs[1,1] + sma2.coefs[1,2]*min(log10_beta_prime[which(spec=='erato')]), yend = sma2.coefs[1,1]+
                     sma2.coefs[1,2]*max(log10_beta_prime[which(spec=='erato')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_beta_prime[which(spec=='iulia')]), xend = max(log10_beta_prime[which(spec=='iulia')]), 
                   y = sma2.coefs[2,1] + sma2.coefs[2,2]*min(log10_beta_prime[which(spec=='iulia')]), yend = sma2.coefs[2,1]+
                     sma2.coefs[2,2]*max(log10_beta_prime[which(spec=='iulia')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_beta_prime[which(spec=='juno')]), xend = max(log10_beta_prime[which(spec=='juno')]), 
                   y = sma2.coefs[3,1] + sma2.coefs[3,2]*min(log10_beta_prime[which(spec=='juno')]), yend = sma2.coefs[3,1]+
                     sma2.coefs[3,2]*max(log10_beta_prime[which(spec=='juno')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_beta_prime[which(spec=='melpomene')]), xend = max(log10_beta_prime[which(spec=='melpomene')]), 
                   y = sma2.coefs[4,1] + sma2.coefs[4,2]*min(log10_beta_prime[which(spec=='melpomene')]), yend = sma2.coefs[4,1]+
                     sma2.coefs[4,2]*max(log10_beta_prime[which(spec=='melpomene')])),colour="darkgrey",linewidth=1)+
  scale_colour_manual(values = c("#c22f20ca", "#1293b9ff", "#7dadb5ff","#cb8781ff"))
lobp2.plot.l

lobp3.plot.l=lobp3.plot+
  geom_segment(aes(x = min(log10_alpha[which(spec=='erato')]), xend = max(log10_alpha[which(spec=='erato')]), 
                   y = sma3.coefs[1,1] + sma3.coefs[1,2]*min(log10_alpha[which(spec=='erato')]), yend = sma3.coefs[1,1]+
                     sma3.coefs[1,2]*max(log10_alpha[which(spec=='erato')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_alpha[which(spec=='iulia')]), xend = max(log10_alpha[which(spec=='iulia')]), 
                   y = sma3.coefs[2,1] + sma3.coefs[2,2]*min(log10_alpha[which(spec=='iulia')]), yend = sma3.coefs[2,1]+
                     sma3.coefs[2,2]*max(log10_alpha[which(spec=='iulia')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_alpha[which(spec=='juno')]), xend = max(log10_alpha[which(spec=='juno')]), 
                   y = sma3.coefs[3,1] + sma3.coefs[3,2]*min(log10_alpha[which(spec=='juno')]), yend = sma3.coefs[3,1]+
                     sma3.coefs[3,2]*max(log10_alpha[which(spec=='juno')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_alpha[which(spec=='melpomene')]), xend = max(log10_alpha[which(spec=='melpomene')]), 
                   y = sma3.coefs[4,1] + sma3.coefs[4,2]*min(log10_alpha[which(spec=='melpomene')]), yend = sma3.coefs[4,1]+
                     sma3.coefs[4,2]*max(log10_alpha[which(spec=='melpomene')])),colour="darkgrey",linewidth=1)+
  scale_colour_manual(values = c("#c22f20ca", "#1293b9ff", "#7dadb5ff","#cb8781ff"))
lobp3.plot.l

lobp4.plot.l=lobp4.plot+
  geom_segment(aes(x = min(log10_beta[which(spec=='erato')]), xend = max(log10_beta[which(spec=='erato')]), 
                   y = sma4.coefs[1,1] + sma4.coefs[1,2]*min(log10_beta[which(spec=='erato')]), yend = sma4.coefs[1,1]+
                     sma4.coefs[1,2]*max(log10_beta[which(spec=='erato')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_beta[which(spec=='iulia')]), xend = max(log10_beta[which(spec=='iulia')]), 
                   y = sma4.coefs[2,1] + sma4.coefs[2,2]*min(log10_beta[which(spec=='iulia')]), yend = sma4.coefs[2,1]+
                     sma4.coefs[2,2]*max(log10_beta[which(spec=='iulia')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_beta[which(spec=='juno')]), xend = max(log10_beta[which(spec=='juno')]), 
                   y = sma4.coefs[3,1] + sma4.coefs[3,2]*min(log10_beta[which(spec=='juno')]), yend = sma4.coefs[3,1]+
                     sma4.coefs[3,2]*max(log10_beta[which(spec=='juno')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_beta[which(spec=='melpomene')]), xend = max(log10_beta[which(spec=='melpomene')]), 
                   y = sma4.coefs[4,1] + sma4.coefs[4,2]*min(log10_beta[which(spec=='melpomene')]), yend = sma4.coefs[4,1]+
                     sma4.coefs[4,2]*max(log10_beta[which(spec=='melpomene')])),colour="darkgrey",linewidth=1)+
  scale_colour_manual(values = c("#c22f20ca", "#1293b9ff", "#7dadb5ff","#cb8781ff"))
lobp4.plot.l

lobp5.plot.l=lobp5.plot+
  geom_segment(aes(x = min(log10_Y.lobe[which(spec=='erato')]), xend = max(log10_Y.lobe[which(spec=='erato')]), 
                   y = sma5.coefs[1,1] + sma5.coefs[1,2]*min(log10_Y.lobe[which(spec=='erato')]), yend = sma5.coefs[1,1]+
                     sma5.coefs[1,2]*max(log10_Y.lobe[which(spec=='erato')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_Y.lobe[which(spec=='iulia')]), xend = max(log10_Y.lobe[which(spec=='iulia')]), 
                   y = sma5.coefs[2,1] + sma5.coefs[2,2]*min(log10_Y.lobe[which(spec=='iulia')]), yend = sma5.coefs[2,1]+
                     sma5.coefs[2,2]*max(log10_Y.lobe[which(spec=='iulia')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_Y.lobe[which(spec=='juno')]), xend = max(log10_Y.lobe[which(spec=='juno')]), 
                   y = sma5.coefs[3,1] + sma5.coefs[3,2]*min(log10_Y.lobe[which(spec=='juno')]), yend = sma5.coefs[3,1]+
                     sma5.coefs[3,2]*max(log10_Y.lobe[which(spec=='juno')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_Y.lobe[which(spec=='melpomene')]), xend = max(log10_Y.lobe[which(spec=='melpomene')]), 
                   y = sma5.coefs[4,1] + sma5.coefs[4,2]*min(log10_Y.lobe[which(spec=='melpomene')]), yend = sma5.coefs[4,1]+
                     sma5.coefs[4,2]*max(log10_Y.lobe[which(spec=='melpomene')])),colour="darkgrey",linewidth=1)+
  scale_colour_manual(values = c("#c22f20ca", "#1293b9ff", "#7dadb5ff","#cb8781ff"))
lobp5.plot.l

lobp6.plot.l=lobp6.plot+
  geom_segment(aes(x = min(log10_gamma[which(spec=='erato')]), xend = max(log10_gamma[which(spec=='erato')]), 
                   y = sma6.coefs[1,1] + sma6.coefs[1,2]*min(log10_gamma[which(spec=='erato')]), yend = sma6.coefs[1,1]+
                     sma6.coefs[1,2]*max(log10_gamma[which(spec=='erato')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_gamma[which(spec=='iulia')]), xend = max(log10_gamma[which(spec=='iulia')]), 
                   y = sma6.coefs[2,1] + sma6.coefs[2,2]*min(log10_gamma[which(spec=='iulia')]), yend = sma6.coefs[2,1]+
                     sma6.coefs[2,2]*max(log10_gamma[which(spec=='iulia')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_gamma[which(spec=='juno')]), xend = max(log10_gamma[which(spec=='juno')]), 
                   y = sma6.coefs[3,1] + sma6.coefs[3,2]*min(log10_gamma[which(spec=='juno')]), yend = sma6.coefs[3,1]+
                     sma6.coefs[3,2]*max(log10_gamma[which(spec=='juno')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_gamma[which(spec=='melpomene')]), xend = max(log10_gamma[which(spec=='melpomene')]), 
                   y = sma6.coefs[4,1] + sma6.coefs[4,2]*min(log10_gamma[which(spec=='melpomene')]), yend = sma6.coefs[4,1]+
                     sma6.coefs[4,2]*max(log10_gamma[which(spec=='melpomene')])),colour="darkgrey",linewidth=1)+
  scale_colour_manual(values = c("#c22f20ca", "#1293b9ff", "#7dadb5ff","#cb8781ff"))
lobp6.plot.l

lobp7.plot.l=lobp7.plot+
  geom_segment(aes(x = min(log10_beta_prime[which(spec=='erato')]), xend = max(log10_beta_prime[which(spec=='erato')]), 
                   y = sma7.coefs[1,1] + sma7.coefs[1,2]*min(log10_beta_prime[which(spec=='erato')]), yend = sma7.coefs[1,1]+
                     sma7.coefs[1,2]*max(log10_beta_prime[which(spec=='erato')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_beta_prime[which(spec=='iulia')]), xend = max(log10_beta_prime[which(spec=='iulia')]), 
                   y = sma7.coefs[2,1] + sma7.coefs[2,2]*min(log10_beta_prime[which(spec=='iulia')]), yend = sma7.coefs[2,1]+
                     sma7.coefs[2,2]*max(log10_beta_prime[which(spec=='iulia')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_beta_prime[which(spec=='juno')]), xend = max(log10_beta_prime[which(spec=='juno')]), 
                   y = sma7.coefs[3,1] + sma7.coefs[3,2]*min(log10_beta_prime[which(spec=='juno')]), yend = sma7.coefs[3,1]+
                     sma7.coefs[3,2]*max(log10_beta_prime[which(spec=='juno')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_beta_prime[which(spec=='melpomene')]), xend = max(log10_beta_prime[which(spec=='melpomene')]), 
                   y = sma7.coefs[4,1] + sma7.coefs[4,2]*min(log10_beta_prime[which(spec=='melpomene')]), yend = sma7.coefs[4,1]+
                     sma7.coefs[4,2]*max(log10_beta_prime[which(spec=='melpomene')])),colour="darkgrey",linewidth=1)+
  scale_colour_manual(values = c("#c22f20ca", "#1293b9ff", "#7dadb5ff","#cb8781ff"))
lobp7.plot.l

lobp8.plot.l=lobp8.plot+
  geom_segment(aes(x = min(log10_alpha[which(spec=='erato')]), xend = max(log10_alpha[which(spec=='erato')]), 
                   y = sma8.coefs[1,1] + sma8.coefs[1,2]*min(log10_alpha[which(spec=='erato')]), yend = sma8.coefs[1,1]+
                     sma8.coefs[1,2]*max(log10_alpha[which(spec=='erato')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_alpha[which(spec=='iulia')]), xend = max(log10_alpha[which(spec=='iulia')]), 
                   y = sma8.coefs[2,1] + sma8.coefs[2,2]*min(log10_alpha[which(spec=='iulia')]), yend = sma8.coefs[2,1]+
                     sma8.coefs[2,2]*max(log10_alpha[which(spec=='iulia')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_alpha[which(spec=='juno')]), xend = max(log10_alpha[which(spec=='juno')]), 
                   y = sma8.coefs[3,1] + sma8.coefs[3,2]*min(log10_alpha[which(spec=='juno')]), yend = sma8.coefs[3,1]+
                     sma8.coefs[3,2]*max(log10_alpha[which(spec=='juno')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_alpha[which(spec=='melpomene')]), xend = max(log10_alpha[which(spec=='melpomene')]), 
                   y = sma8.coefs[4,1] + sma8.coefs[4,2]*min(log10_alpha[which(spec=='melpomene')]), yend = sma8.coefs[4,1]+
                     sma8.coefs[4,2]*max(log10_alpha[which(spec=='melpomene')])),colour="darkgrey",linewidth=1)+
  scale_colour_manual(values = c("#c22f20ca", "#1293b9ff", "#7dadb5ff","#cb8781ff"))
lobp8.plot.l

lobp9.plot.l=lobp9.plot+
  geom_segment(aes(x = min(log10_beta[which(spec=='erato')]), xend = max(log10_beta[which(spec=='erato')]), 
                   y = sma9.coefs[1,1] + sma9.coefs[1,2]*min(log10_beta[which(spec=='erato')]), yend = sma9.coefs[1,1]+
                     sma9.coefs[1,2]*max(log10_beta[which(spec=='erato')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_beta[which(spec=='iulia')]), xend = max(log10_beta[which(spec=='iulia')]), 
                   y = sma9.coefs[2,1] + sma9.coefs[2,2]*min(log10_beta[which(spec=='iulia')]), yend = sma9.coefs[2,1]+
                     sma9.coefs[2,2]*max(log10_beta[which(spec=='iulia')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_beta[which(spec=='juno')]), xend = max(log10_beta[which(spec=='juno')]), 
                   y = sma9.coefs[3,1] + sma9.coefs[3,2]*min(log10_beta[which(spec=='juno')]), yend = sma9.coefs[3,1]+
                     sma9.coefs[3,2]*max(log10_beta[which(spec=='juno')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_beta[which(spec=='melpomene')]), xend = max(log10_beta[which(spec=='melpomene')]), 
                   y = sma9.coefs[4,1] + sma9.coefs[4,2]*min(log10_beta[which(spec=='melpomene')]), yend = sma9.coefs[4,1]+
                     sma9.coefs[4,2]*max(log10_beta[which(spec=='melpomene')])),colour="darkgrey",linewidth=1)+
  scale_colour_manual(values = c("#c22f20ca", "#1293b9ff", "#7dadb5ff","#cb8781ff"))
lobp9.plot.l

lobp10.plot.l=lobp10.plot+
  geom_segment(aes(x = min(log10_Y.lobe[which(spec=='erato')]), xend = max(log10_Y.lobe[which(spec=='erato')]), 
                   y = sma10.coefs[1,1] + sma10.coefs[1,2]*min(log10_Y.lobe[which(spec=='erato')]), yend = sma10.coefs[1,1]+
                     sma10.coefs[1,2]*max(log10_Y.lobe[which(spec=='erato')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_Y.lobe[which(spec=='iulia')]), xend = max(log10_Y.lobe[which(spec=='iulia')]), 
                   y = sma10.coefs[2,1] + sma10.coefs[2,2]*min(log10_Y.lobe[which(spec=='iulia')]), yend = sma10.coefs[2,1]+
                     sma10.coefs[2,2]*max(log10_Y.lobe[which(spec=='iulia')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_Y.lobe[which(spec=='juno')]), xend = max(log10_Y.lobe[which(spec=='juno')]), 
                   y = sma10.coefs[3,1] + sma10.coefs[3,2]*min(log10_Y.lobe[which(spec=='juno')]), yend = sma10.coefs[3,1]+
                     sma10.coefs[3,2]*max(log10_Y.lobe[which(spec=='juno')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_Y.lobe[which(spec=='melpomene')]), xend = max(log10_Y.lobe[which(spec=='melpomene')]), 
                   y = sma10.coefs[4,1] + sma10.coefs[4,2]*min(log10_Y.lobe[which(spec=='melpomene')]), yend = sma10.coefs[4,1]+
                     sma10.coefs[4,2]*max(log10_Y.lobe[which(spec=='melpomene')])),colour="darkgrey",linewidth=1)+
  scale_colour_manual(values = c("#c22f20ca", "#1293b9ff", "#7dadb5ff","#cb8781ff"))
lobp10.plot.l

lobp11.plot.l=lobp11.plot+
  geom_segment(aes(x = min(log10_gamma[which(spec=='erato')]), xend = max(log10_gamma[which(spec=='erato')]), 
                   y = sma11.coefs[1,1] + sma11.coefs[1,2]*min(log10_gamma[which(spec=='erato')]), yend = sma11.coefs[1,1]+
                     sma11.coefs[1,2]*max(log10_gamma[which(spec=='erato')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_gamma[which(spec=='iulia')]), xend = max(log10_gamma[which(spec=='iulia')]), 
                   y = sma11.coefs[2,1] + sma11.coefs[2,2]*min(log10_gamma[which(spec=='iulia')]), yend = sma11.coefs[2,1]+
                     sma11.coefs[2,2]*max(log10_gamma[which(spec=='iulia')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_gamma[which(spec=='juno')]), xend = max(log10_gamma[which(spec=='juno')]), 
                   y = sma11.coefs[3,1] + sma11.coefs[3,2]*min(log10_gamma[which(spec=='juno')]), yend = sma11.coefs[3,1]+
                     sma11.coefs[3,2]*max(log10_gamma[which(spec=='juno')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_gamma[which(spec=='melpomene')]), xend = max(log10_gamma[which(spec=='melpomene')]), 
                   y = sma11.coefs[4,1] + sma11.coefs[4,2]*min(log10_gamma[which(spec=='melpomene')]), yend = sma11.coefs[4,1]+
                     sma11.coefs[4,2]*max(log10_gamma[which(spec=='melpomene')])),colour="darkgrey",linewidth=1)+
  scale_colour_manual(values = c("#c22f20ca", "#1293b9ff", "#7dadb5ff","#cb8781ff"))
lobp11.plot.l

lobp12.plot.l=lobp12.plot+
  geom_segment(aes(x = min(log10_alpha_prime[which(spec=='erato')]), xend = max(log10_alpha_prime[which(spec=='erato')]), 
                   y = sma12.coefs[1,1] + sma12.coefs[1,2]*min(log10_alpha_prime[which(spec=='erato')]), yend = sma12.coefs[1,1]+
                     sma12.coefs[1,2]*max(log10_alpha_prime[which(spec=='erato')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_alpha_prime[which(spec=='iulia')]), xend = max(log10_alpha_prime[which(spec=='iulia')]), 
                   y = sma12.coefs[2,1] + sma12.coefs[2,2]*min(log10_alpha_prime[which(spec=='iulia')]), yend = sma12.coefs[2,1]+
                     sma12.coefs[2,2]*max(log10_alpha_prime[which(spec=='iulia')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_alpha_prime[which(spec=='juno')]), xend = max(log10_alpha_prime[which(spec=='juno')]), 
                   y = sma12.coefs[3,1] + sma12.coefs[3,2]*min(log10_alpha_prime[which(spec=='juno')]), yend = sma12.coefs[3,1]+
                     sma12.coefs[3,2]*max(log10_alpha_prime[which(spec=='juno')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_alpha_prime[which(spec=='melpomene')]), xend = max(log10_alpha_prime[which(spec=='melpomene')]), 
                   y = sma12.coefs[4,1] + sma12.coefs[4,2]*min(log10_alpha_prime[which(spec=='melpomene')]), yend = sma12.coefs[4,1]+
                     sma12.coefs[4,2]*max(log10_alpha_prime[which(spec=='melpomene')])),colour="darkgrey",linewidth=1)+
  scale_colour_manual(values = c("#c22f20ca", "#1293b9ff", "#7dadb5ff","#cb8781ff"))
lobp12.plot.l

lobp13.plot.l=lobp13.plot+
  geom_segment(aes(x = min(log10_alpha[which(spec=='erato')]), xend = max(log10_alpha[which(spec=='erato')]), 
                   y = sma13.coefs[1,1] + sma13.coefs[1,2]*min(log10_alpha[which(spec=='erato')]), yend = sma13.coefs[1,1]+
                     sma13.coefs[1,2]*max(log10_alpha[which(spec=='erato')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_alpha[which(spec=='iulia')]), xend = max(log10_alpha[which(spec=='iulia')]), 
                   y = sma13.coefs[2,1] + sma13.coefs[2,2]*min(log10_alpha[which(spec=='iulia')]), yend = sma13.coefs[2,1]+
                     sma13.coefs[2,2]*max(log10_alpha[which(spec=='iulia')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_alpha[which(spec=='juno')]), xend = max(log10_alpha[which(spec=='juno')]), 
                   y = sma13.coefs[3,1] + sma13.coefs[3,2]*min(log10_alpha[which(spec=='juno')]), yend = sma13.coefs[3,1]+
                     sma13.coefs[3,2]*max(log10_alpha[which(spec=='juno')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_alpha[which(spec=='melpomene')]), xend = max(log10_alpha[which(spec=='melpomene')]), 
                   y = sma13.coefs[4,1] + sma13.coefs[4,2]*min(log10_alpha[which(spec=='melpomene')]), yend = sma13.coefs[4,1]+
                     sma13.coefs[4,2]*max(log10_alpha[which(spec=='melpomene')])),colour="darkgrey",linewidth=1)+
  scale_colour_manual(values = c("#c22f20ca", "#1293b9ff", "#7dadb5ff","#cb8781ff"))
lobp13.plot.l

lobp14.plot.l=lobp14.plot+
  geom_segment(aes(x = min(log10_beta[which(spec=='erato')]), xend = max(log10_beta[which(spec=='erato')]), 
                   y = sma14.coefs[1,1] + sma14.coefs[1,2]*min(log10_beta[which(spec=='erato')]), yend = sma14.coefs[1,1]+
                     sma14.coefs[1,2]*max(log10_beta[which(spec=='erato')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_beta[which(spec=='iulia')]), xend = max(log10_beta[which(spec=='iulia')]), 
                   y = sma14.coefs[2,1] + sma14.coefs[2,2]*min(log10_beta[which(spec=='iulia')]), yend = sma14.coefs[2,1]+
                     sma14.coefs[2,2]*max(log10_beta[which(spec=='iulia')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_beta[which(spec=='juno')]), xend = max(log10_beta[which(spec=='juno')]), 
                   y = sma14.coefs[3,1] + sma14.coefs[3,2]*min(log10_beta[which(spec=='juno')]), yend = sma14.coefs[3,1]+
                     sma14.coefs[3,2]*max(log10_beta[which(spec=='juno')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_beta[which(spec=='melpomene')]), xend = max(log10_beta[which(spec=='melpomene')]), 
                   y = sma14.coefs[4,1] + sma14.coefs[4,2]*min(log10_beta[which(spec=='melpomene')]), yend = sma14.coefs[4,1]+
                     sma14.coefs[4,2]*max(log10_beta[which(spec=='melpomene')])),colour="darkgrey",linewidth=1)+
  scale_colour_manual(values = c("#c22f20ca", "#1293b9ff", "#7dadb5ff","#cb8781ff"))
lobp14.plot.l

lobp15.plot.l=lobp15.plot+
  geom_segment(aes(x = min(log10_Y.lobe[which(spec=='erato')]), xend = max(log10_Y.lobe[which(spec=='erato')]), 
                   y = sma15.coefs[1,1] + sma15.coefs[1,2]*min(log10_Y.lobe[which(spec=='erato')]), yend = sma15.coefs[1,1]+
                     sma15.coefs[1,2]*max(log10_Y.lobe[which(spec=='erato')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_Y.lobe[which(spec=='iulia')]), xend = max(log10_Y.lobe[which(spec=='iulia')]), 
                   y = sma15.coefs[2,1] + sma15.coefs[2,2]*min(log10_Y.lobe[which(spec=='iulia')]), yend = sma15.coefs[2,1]+
                     sma15.coefs[2,2]*max(log10_Y.lobe[which(spec=='iulia')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_Y.lobe[which(spec=='juno')]), xend = max(log10_Y.lobe[which(spec=='juno')]), 
                   y = sma15.coefs[3,1] + sma15.coefs[3,2]*min(log10_Y.lobe[which(spec=='juno')]), yend = sma15.coefs[3,1]+
                     sma15.coefs[3,2]*max(log10_Y.lobe[which(spec=='juno')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_Y.lobe[which(spec=='melpomene')]), xend = max(log10_Y.lobe[which(spec=='melpomene')]), 
                   y = sma15.coefs[4,1] + sma15.coefs[4,2]*min(log10_Y.lobe[which(spec=='melpomene')]), yend = sma15.coefs[4,1]+
                     sma15.coefs[4,2]*max(log10_Y.lobe[which(spec=='melpomene')])),colour="darkgrey",linewidth=1)+
  scale_colour_manual(values = c("#c22f20ca", "#1293b9ff", "#7dadb5ff","#cb8781ff"))
lobp15.plot.l

lobp16.plot.l=lobp16.plot+
  geom_segment(aes(x = min(log10_gamma[which(spec=='erato')]), xend = max(log10_gamma[which(spec=='erato')]), 
                   y = sma16.coefs[1,1] + sma16.coefs[1,2]*min(log10_gamma[which(spec=='erato')]), yend = sma16.coefs[1,1]+
                     sma16.coefs[1,2]*max(log10_gamma[which(spec=='erato')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_gamma[which(spec=='iulia')]), xend = max(log10_gamma[which(spec=='iulia')]), 
                   y = sma16.coefs[2,1] + sma16.coefs[2,2]*min(log10_gamma[which(spec=='iulia')]), yend = sma16.coefs[2,1]+
                     sma16.coefs[2,2]*max(log10_gamma[which(spec=='iulia')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_gamma[which(spec=='juno')]), xend = max(log10_gamma[which(spec=='juno')]), 
                   y = sma16.coefs[3,1] + sma16.coefs[3,2]*min(log10_gamma[which(spec=='juno')]), yend = sma16.coefs[3,1]+
                     sma16.coefs[3,2]*max(log10_gamma[which(spec=='juno')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_gamma[which(spec=='melpomene')]), xend = max(log10_gamma[which(spec=='melpomene')]), 
                   y = sma16.coefs[4,1] + sma16.coefs[4,2]*min(log10_gamma[which(spec=='melpomene')]), yend = sma16.coefs[4,1]+
                     sma16.coefs[4,2]*max(log10_gamma[which(spec=='melpomene')])),colour="darkgrey",linewidth=1)+
  scale_colour_manual(values = c("#c22f20ca", "#1293b9ff", "#7dadb5ff","#cb8781ff"))
lobp16.plot.l

lobp17.plot.l=lobp17.plot+
  geom_segment(aes(x = min(log10_alpha_prime[which(spec=='erato')]), xend = max(log10_alpha_prime[which(spec=='erato')]), 
                   y = sma17.coefs[1,1] + sma17.coefs[1,2]*min(log10_alpha_prime[which(spec=='erato')]), yend = sma17.coefs[1,1]+
                     sma17.coefs[1,2]*max(log10_alpha_prime[which(spec=='erato')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_alpha_prime[which(spec=='iulia')]), xend = max(log10_alpha_prime[which(spec=='iulia')]), 
                   y = sma17.coefs[2,1] + sma17.coefs[2,2]*min(log10_alpha_prime[which(spec=='iulia')]), yend = sma17.coefs[2,1]+
                     sma17.coefs[2,2]*max(log10_alpha_prime[which(spec=='iulia')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_alpha_prime[which(spec=='juno')]), xend = max(log10_alpha_prime[which(spec=='juno')]), 
                   y = sma17.coefs[3,1] + sma17.coefs[3,2]*min(log10_alpha_prime[which(spec=='juno')]), yend = sma17.coefs[3,1]+
                     sma17.coefs[3,2]*max(log10_alpha_prime[which(spec=='juno')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_alpha_prime[which(spec=='melpomene')]), xend = max(log10_alpha_prime[which(spec=='melpomene')]), 
                   y = sma17.coefs[4,1] + sma17.coefs[4,2]*min(log10_alpha_prime[which(spec=='melpomene')]), yend = sma17.coefs[4,1]+
                     sma17.coefs[4,2]*max(log10_alpha_prime[which(spec=='melpomene')])),colour="darkgrey",linewidth=1)+
  scale_colour_manual(values = c("#c22f20ca", "#1293b9ff", "#7dadb5ff","#cb8781ff"))
lobp17.plot.l

lobp18.plot.l=lobp18.plot+
  geom_segment(aes(x = min(log10_beta_prime[which(spec=='erato')]), xend = max(log10_beta_prime[which(spec=='erato')]), 
                   y = sma18.coefs[1,1] + sma18.coefs[1,2]*min(log10_beta_prime[which(spec=='erato')]), yend = sma18.coefs[1,1]+
                     sma18.coefs[1,2]*max(log10_beta_prime[which(spec=='erato')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_beta_prime[which(spec=='iulia')]), xend = max(log10_beta_prime[which(spec=='iulia')]), 
                   y = sma18.coefs[2,1] + sma18.coefs[2,2]*min(log10_beta_prime[which(spec=='iulia')]), yend = sma18.coefs[2,1]+
                     sma18.coefs[2,2]*max(log10_beta_prime[which(spec=='iulia')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_beta_prime[which(spec=='juno')]), xend = max(log10_beta_prime[which(spec=='juno')]), 
                   y = sma18.coefs[3,1] + sma18.coefs[3,2]*min(log10_beta_prime[which(spec=='juno')]), yend = sma18.coefs[3,1]+
                     sma18.coefs[3,2]*max(log10_beta_prime[which(spec=='juno')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_beta_prime[which(spec=='melpomene')]), xend = max(log10_beta_prime[which(spec=='melpomene')]), 
                   y = sma18.coefs[4,1] + sma18.coefs[4,2]*min(log10_beta_prime[which(spec=='melpomene')]), yend = sma18.coefs[4,1]+
                     sma18.coefs[4,2]*max(log10_beta_prime[which(spec=='melpomene')])),colour="darkgrey",linewidth=1)+
  scale_colour_manual(values = c("#c22f20ca", "#1293b9ff", "#7dadb5ff","#cb8781ff"))
lobp18.plot.l

lobp19.plot.l=lobp19.plot+
  geom_segment(aes(x = min(log10_beta[which(spec=='erato')]), xend = max(log10_beta[which(spec=='erato')]), 
                   y = sma19.coefs[1,1] + sma19.coefs[1,2]*min(log10_beta[which(spec=='erato')]), yend = sma19.coefs[1,1]+
                     sma19.coefs[1,2]*max(log10_beta[which(spec=='erato')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_beta[which(spec=='iulia')]), xend = max(log10_beta[which(spec=='iulia')]), 
                   y = sma19.coefs[2,1] + sma19.coefs[2,2]*min(log10_beta[which(spec=='iulia')]), yend = sma19.coefs[2,1]+
                     sma19.coefs[2,2]*max(log10_beta[which(spec=='iulia')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_beta[which(spec=='juno')]), xend = max(log10_beta[which(spec=='juno')]), 
                   y = sma19.coefs[3,1] + sma19.coefs[3,2]*min(log10_beta[which(spec=='juno')]), yend = sma19.coefs[3,1]+
                     sma19.coefs[3,2]*max(log10_beta[which(spec=='juno')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_beta[which(spec=='melpomene')]), xend = max(log10_beta[which(spec=='melpomene')]), 
                   y = sma19.coefs[4,1] + sma19.coefs[4,2]*min(log10_beta[which(spec=='melpomene')]), yend = sma19.coefs[4,1]+
                     sma19.coefs[4,2]*max(log10_beta[which(spec=='melpomene')])),colour="darkgrey",linewidth=1)+
  scale_colour_manual(values = c("#c22f20ca", "#1293b9ff", "#7dadb5ff","#cb8781ff"))
lobp19.plot.l

lobp20.plot.l=lobp20.plot+
  geom_segment(aes(x = min(log10_Y.lobe[which(spec=='erato')]), xend = max(log10_Y.lobe[which(spec=='erato')]), 
                   y = sma20.coefs[1,1] + sma20.coefs[1,2]*min(log10_Y.lobe[which(spec=='erato')]), yend = sma20.coefs[1,1]+
                     sma20.coefs[1,2]*max(log10_Y.lobe[which(spec=='erato')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_Y.lobe[which(spec=='iulia')]), xend = max(log10_Y.lobe[which(spec=='iulia')]), 
                   y = sma20.coefs[2,1] + sma20.coefs[2,2]*min(log10_Y.lobe[which(spec=='iulia')]), yend = sma20.coefs[2,1]+
                     sma20.coefs[2,2]*max(log10_Y.lobe[which(spec=='iulia')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_Y.lobe[which(spec=='juno')]), xend = max(log10_Y.lobe[which(spec=='juno')]), 
                   y = sma20.coefs[3,1] + sma20.coefs[3,2]*min(log10_Y.lobe[which(spec=='juno')]), yend = sma20.coefs[3,1]+
                     sma20.coefs[3,2]*max(log10_Y.lobe[which(spec=='juno')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_Y.lobe[which(spec=='melpomene')]), xend = max(log10_Y.lobe[which(spec=='melpomene')]), 
                   y = sma20.coefs[4,1] + sma20.coefs[4,2]*min(log10_Y.lobe[which(spec=='melpomene')]), yend = sma20.coefs[4,1]+
                     sma20.coefs[4,2]*max(log10_Y.lobe[which(spec=='melpomene')])),colour="darkgrey",linewidth=1)+
  scale_colour_manual(values = c("#c22f20ca", "#1293b9ff", "#7dadb5ff","#cb8781ff"))
lobp20.plot.l

lobp21.plot.l=lobp21.plot+
  geom_segment(aes(x = min(log10_gamma[which(spec=='erato')]), xend = max(log10_gamma[which(spec=='erato')]), 
                   y = sma21.coefs[1,1] + sma21.coefs[1,2]*min(log10_gamma[which(spec=='erato')]), yend = sma21.coefs[1,1]+
                     sma21.coefs[1,2]*max(log10_gamma[which(spec=='erato')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_gamma[which(spec=='iulia')]), xend = max(log10_gamma[which(spec=='iulia')]), 
                   y = sma21.coefs[2,1] + sma21.coefs[2,2]*min(log10_gamma[which(spec=='iulia')]), yend = sma21.coefs[2,1]+
                     sma21.coefs[2,2]*max(log10_gamma[which(spec=='iulia')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_gamma[which(spec=='juno')]), xend = max(log10_gamma[which(spec=='juno')]), 
                   y = sma21.coefs[3,1] + sma21.coefs[3,2]*min(log10_gamma[which(spec=='juno')]), yend = sma21.coefs[3,1]+
                     sma21.coefs[3,2]*max(log10_gamma[which(spec=='juno')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_gamma[which(spec=='melpomene')]), xend = max(log10_gamma[which(spec=='melpomene')]), 
                   y = sma21.coefs[4,1] + sma21.coefs[4,2]*min(log10_gamma[which(spec=='melpomene')]), yend = sma21.coefs[4,1]+
                     sma21.coefs[4,2]*max(log10_gamma[which(spec=='melpomene')])),colour="darkgrey",linewidth=1)+
  scale_colour_manual(values = c("#c22f20ca", "#1293b9ff", "#7dadb5ff","#cb8781ff"))
lobp21.plot.l

lobp22.plot.l=lobp22.plot+
  geom_segment(aes(x = min(log10_alpha_prime[which(spec=='erato')]), xend = max(log10_alpha_prime[which(spec=='erato')]), 
                   y = sma22.coefs[1,1] + sma22.coefs[1,2]*min(log10_alpha_prime[which(spec=='erato')]), yend = sma22.coefs[1,1]+
                     sma22.coefs[1,2]*max(log10_alpha_prime[which(spec=='erato')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_alpha_prime[which(spec=='iulia')]), xend = max(log10_alpha_prime[which(spec=='iulia')]), 
                   y = sma22.coefs[2,1] + sma22.coefs[2,2]*min(log10_alpha_prime[which(spec=='iulia')]), yend = sma22.coefs[2,1]+
                     sma22.coefs[2,2]*max(log10_alpha_prime[which(spec=='iulia')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_alpha_prime[which(spec=='juno')]), xend = max(log10_alpha_prime[which(spec=='juno')]), 
                   y = sma22.coefs[3,1] + sma22.coefs[3,2]*min(log10_alpha_prime[which(spec=='juno')]), yend = sma22.coefs[3,1]+
                     sma22.coefs[3,2]*max(log10_alpha_prime[which(spec=='juno')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_alpha_prime[which(spec=='melpomene')]), xend = max(log10_alpha_prime[which(spec=='melpomene')]), 
                   y = sma22.coefs[4,1] + sma22.coefs[4,2]*min(log10_alpha_prime[which(spec=='melpomene')]), yend = sma22.coefs[4,1]+
                     sma22.coefs[4,2]*max(log10_alpha_prime[which(spec=='melpomene')])),colour="darkgrey",linewidth=1)+
  scale_colour_manual(values = c("#c22f20ca", "#1293b9ff", "#7dadb5ff","#cb8781ff"))
lobp22.plot.l

lobp23.plot.l=lobp23.plot+
  geom_segment(aes(x = min(log10_beta_prime[which(spec=='erato')]), xend = max(log10_beta_prime[which(spec=='erato')]), 
                   y = sma23.coefs[1,1] + sma23.coefs[1,2]*min(log10_beta_prime[which(spec=='erato')]), yend = sma23.coefs[1,1]+
                     sma23.coefs[1,2]*max(log10_beta_prime[which(spec=='erato')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_beta_prime[which(spec=='iulia')]), xend = max(log10_beta_prime[which(spec=='iulia')]), 
                   y = sma23.coefs[2,1] + sma23.coefs[2,2]*min(log10_beta_prime[which(spec=='iulia')]), yend = sma23.coefs[2,1]+
                     sma23.coefs[2,2]*max(log10_beta_prime[which(spec=='iulia')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_beta_prime[which(spec=='juno')]), xend = max(log10_beta_prime[which(spec=='juno')]), 
                   y = sma23.coefs[3,1] + sma23.coefs[3,2]*min(log10_beta_prime[which(spec=='juno')]), yend = sma23.coefs[3,1]+
                     sma23.coefs[3,2]*max(log10_beta_prime[which(spec=='juno')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_beta_prime[which(spec=='melpomene')]), xend = max(log10_beta_prime[which(spec=='melpomene')]), 
                   y = sma23.coefs[4,1] + sma23.coefs[4,2]*min(log10_beta_prime[which(spec=='melpomene')]), yend = sma23.coefs[4,1]+
                     sma23.coefs[4,2]*max(log10_beta_prime[which(spec=='melpomene')])),colour="darkgrey",linewidth=1)+
  scale_colour_manual(values = c("#c22f20ca", "#1293b9ff", "#7dadb5ff","#cb8781ff"))
lobp23.plot.l

lobp24.plot.l=lobp24.plot+
  geom_segment(aes(x = min(log10_alpha[which(spec=='erato')]), xend = max(log10_alpha[which(spec=='erato')]), 
                   y = sma24.coefs[1,1] + sma24.coefs[1,2]*min(log10_alpha[which(spec=='erato')]), yend = sma24.coefs[1,1]+
                     sma24.coefs[1,2]*max(log10_alpha[which(spec=='erato')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_alpha[which(spec=='iulia')]), xend = max(log10_alpha[which(spec=='iulia')]), 
                   y = sma24.coefs[2,1] + sma24.coefs[2,2]*min(log10_alpha[which(spec=='iulia')]), yend = sma24.coefs[2,1]+
                     sma24.coefs[2,2]*max(log10_alpha[which(spec=='iulia')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_alpha[which(spec=='juno')]), xend = max(log10_alpha[which(spec=='juno')]), 
                   y = sma24.coefs[3,1] + sma24.coefs[3,2]*min(log10_alpha[which(spec=='juno')]), yend = sma24.coefs[3,1]+
                     sma24.coefs[3,2]*max(log10_alpha[which(spec=='juno')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_alpha[which(spec=='melpomene')]), xend = max(log10_alpha[which(spec=='melpomene')]), 
                   y = sma24.coefs[4,1] + sma24.coefs[4,2]*min(log10_alpha[which(spec=='melpomene')]), yend = sma24.coefs[4,1]+
                     sma24.coefs[4,2]*max(log10_alpha[which(spec=='melpomene')])),colour="darkgrey",linewidth=1)+
  scale_colour_manual(values = c("#c22f20ca", "#1293b9ff", "#7dadb5ff","#cb8781ff"))
lobp24.plot.l

lobp25.plot.l=lobp25.plot+
  geom_segment(aes(x = min(log10_Y.lobe[which(spec=='erato')]), xend = max(log10_Y.lobe[which(spec=='erato')]), 
                   y = sma25.coefs[1,1] + sma25.coefs[1,2]*min(log10_Y.lobe[which(spec=='erato')]), yend = sma25.coefs[1,1]+
                     sma25.coefs[1,2]*max(log10_Y.lobe[which(spec=='erato')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_Y.lobe[which(spec=='iulia')]), xend = max(log10_Y.lobe[which(spec=='iulia')]), 
                   y = sma25.coefs[2,1] + sma25.coefs[2,2]*min(log10_Y.lobe[which(spec=='iulia')]), yend = sma25.coefs[2,1]+
                     sma25.coefs[2,2]*max(log10_Y.lobe[which(spec=='iulia')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_Y.lobe[which(spec=='juno')]), xend = max(log10_Y.lobe[which(spec=='juno')]), 
                   y = sma25.coefs[3,1] + sma25.coefs[3,2]*min(log10_Y.lobe[which(spec=='juno')]), yend = sma25.coefs[3,1]+
                     sma25.coefs[3,2]*max(log10_Y.lobe[which(spec=='juno')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_Y.lobe[which(spec=='melpomene')]), xend = max(log10_Y.lobe[which(spec=='melpomene')]), 
                   y = sma25.coefs[4,1] + sma25.coefs[4,2]*min(log10_Y.lobe[which(spec=='melpomene')]), yend = sma25.coefs[4,1]+
                     sma25.coefs[4,2]*max(log10_Y.lobe[which(spec=='melpomene')])),colour="darkgrey",linewidth=1)+
  scale_colour_manual(values = c("#c22f20ca", "#1293b9ff", "#7dadb5ff","#cb8781ff"))
lobp25.plot.l

lobp26.plot.l=lobp26.plot+
  geom_segment(aes(x = min(log10_gamma[which(spec=='erato')]), xend = max(log10_gamma[which(spec=='erato')]), 
                   y = sma26.coefs[1,1] + sma26.coefs[1,2]*min(log10_gamma[which(spec=='erato')]), yend = sma26.coefs[1,1]+
                     sma26.coefs[1,2]*max(log10_gamma[which(spec=='erato')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_gamma[which(spec=='iulia')]), xend = max(log10_gamma[which(spec=='iulia')]), 
                   y = sma26.coefs[2,1] + sma26.coefs[2,2]*min(log10_gamma[which(spec=='iulia')]), yend = sma26.coefs[2,1]+
                     sma26.coefs[2,2]*max(log10_gamma[which(spec=='iulia')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_gamma[which(spec=='juno')]), xend = max(log10_gamma[which(spec=='juno')]), 
                   y = sma26.coefs[3,1] + sma26.coefs[3,2]*min(log10_gamma[which(spec=='juno')]), yend = sma26.coefs[3,1]+
                     sma26.coefs[3,2]*max(log10_gamma[which(spec=='juno')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_gamma[which(spec=='melpomene')]), xend = max(log10_gamma[which(spec=='melpomene')]), 
                   y = sma26.coefs[4,1] + sma26.coefs[4,2]*min(log10_gamma[which(spec=='melpomene')]), yend = sma26.coefs[4,1]+
                     sma26.coefs[4,2]*max(log10_gamma[which(spec=='melpomene')])),colour="darkgrey",linewidth=1)+
  scale_colour_manual(values = c("#c22f20ca", "#1293b9ff", "#7dadb5ff","#cb8781ff"))
lobp26.plot.l

lobp27.plot.l=lobp27.plot+
  geom_segment(aes(x = min(log10_alpha_prime[which(spec=='erato')]), xend = max(log10_alpha_prime[which(spec=='erato')]), 
                   y = sma27.coefs[1,1] + sma27.coefs[1,2]*min(log10_alpha_prime[which(spec=='erato')]), yend = sma27.coefs[1,1]+
                     sma27.coefs[1,2]*max(log10_alpha_prime[which(spec=='erato')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_alpha_prime[which(spec=='iulia')]), xend = max(log10_alpha_prime[which(spec=='iulia')]), 
                   y = sma27.coefs[2,1] + sma27.coefs[2,2]*min(log10_alpha_prime[which(spec=='iulia')]), yend = sma27.coefs[2,1]+
                     sma27.coefs[2,2]*max(log10_alpha_prime[which(spec=='iulia')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_alpha_prime[which(spec=='juno')]), xend = max(log10_alpha_prime[which(spec=='juno')]), 
                   y = sma27.coefs[3,1] + sma27.coefs[3,2]*min(log10_alpha_prime[which(spec=='juno')]), yend = sma27.coefs[3,1]+
                     sma27.coefs[3,2]*max(log10_alpha_prime[which(spec=='juno')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_alpha_prime[which(spec=='melpomene')]), xend = max(log10_alpha_prime[which(spec=='melpomene')]), 
                   y = sma27.coefs[4,1] + sma27.coefs[4,2]*min(log10_alpha_prime[which(spec=='melpomene')]), yend = sma27.coefs[4,1]+
                     sma27.coefs[4,2]*max(log10_alpha_prime[which(spec=='melpomene')])),colour="darkgrey",linewidth=1)+
  scale_colour_manual(values = c("#c22f20ca", "#1293b9ff", "#7dadb5ff","#cb8781ff"))
lobp27.plot.l

lobp28.plot.l=lobp28.plot+
  geom_segment(aes(x = min(log10_beta_prime[which(spec=='erato')]), xend = max(log10_beta_prime[which(spec=='erato')]), 
                   y = sma28.coefs[1,1] + sma28.coefs[1,2]*min(log10_beta_prime[which(spec=='erato')]), yend = sma28.coefs[1,1]+
                     sma28.coefs[1,2]*max(log10_beta_prime[which(spec=='erato')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_beta_prime[which(spec=='iulia')]), xend = max(log10_beta_prime[which(spec=='iulia')]), 
                   y = sma28.coefs[2,1] + sma28.coefs[2,2]*min(log10_beta_prime[which(spec=='iulia')]), yend = sma28.coefs[2,1]+
                     sma28.coefs[2,2]*max(log10_beta_prime[which(spec=='iulia')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_beta_prime[which(spec=='juno')]), xend = max(log10_beta_prime[which(spec=='juno')]), 
                   y = sma28.coefs[3,1] + sma28.coefs[3,2]*min(log10_beta_prime[which(spec=='juno')]), yend = sma28.coefs[3,1]+
                     sma28.coefs[3,2]*max(log10_beta_prime[which(spec=='juno')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_beta_prime[which(spec=='melpomene')]), xend = max(log10_beta_prime[which(spec=='melpomene')]), 
                   y = sma28.coefs[4,1] + sma28.coefs[4,2]*min(log10_beta_prime[which(spec=='melpomene')]), yend = sma28.coefs[4,1]+
                     sma28.coefs[4,2]*max(log10_beta_prime[which(spec=='melpomene')])),colour="darkgrey",linewidth=1)+
  scale_colour_manual(values = c("#c22f20ca", "#1293b9ff", "#7dadb5ff","#cb8781ff"))
lobp28.plot.l

lobp29.plot.l=lobp29.plot+
  geom_segment(aes(x = min(log10_alpha[which(spec=='erato')]), xend = max(log10_alpha[which(spec=='erato')]), 
                   y = sma29.coefs[1,1] + sma29.coefs[1,2]*min(log10_alpha[which(spec=='erato')]), yend = sma29.coefs[1,1]+
                     sma29.coefs[1,2]*max(log10_alpha[which(spec=='erato')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_alpha[which(spec=='iulia')]), xend = max(log10_alpha[which(spec=='iulia')]), 
                   y = sma29.coefs[2,1] + sma29.coefs[2,2]*min(log10_alpha[which(spec=='iulia')]), yend = sma29.coefs[2,1]+
                     sma29.coefs[2,2]*max(log10_alpha[which(spec=='iulia')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_alpha[which(spec=='juno')]), xend = max(log10_alpha[which(spec=='juno')]), 
                   y = sma29.coefs[3,1] + sma29.coefs[3,2]*min(log10_alpha[which(spec=='juno')]), yend = sma29.coefs[3,1]+
                     sma29.coefs[3,2]*max(log10_alpha[which(spec=='juno')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_alpha[which(spec=='melpomene')]), xend = max(log10_alpha[which(spec=='melpomene')]), 
                   y = sma29.coefs[4,1] + sma29.coefs[4,2]*min(log10_alpha[which(spec=='melpomene')]), yend = sma29.coefs[4,1]+
                     sma29.coefs[4,2]*max(log10_alpha[which(spec=='melpomene')])),colour="darkgrey",linewidth=1)+
  scale_colour_manual(values = c("#c22f20ca", "#1293b9ff", "#7dadb5ff","#cb8781ff"))
lobp29.plot.l

lobp30.plot.l=lobp30.plot+
  geom_segment(aes(x = min(log10_beta[which(spec=='erato')]), xend = max(log10_beta[which(spec=='erato')]), 
                   y = sma30.coefs[1,1] + sma30.coefs[1,2]*min(log10_beta[which(spec=='erato')]), yend = sma30.coefs[1,1]+
                     sma30.coefs[1,2]*max(log10_beta[which(spec=='erato')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_beta[which(spec=='iulia')]), xend = max(log10_beta[which(spec=='iulia')]), 
                   y = sma30.coefs[2,1] + sma30.coefs[2,2]*min(log10_beta[which(spec=='iulia')]), yend = sma30.coefs[2,1]+
                     sma30.coefs[2,2]*max(log10_beta[which(spec=='iulia')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_beta[which(spec=='juno')]), xend = max(log10_beta[which(spec=='juno')]), 
                   y = sma30.coefs[3,1] + sma30.coefs[3,2]*min(log10_beta[which(spec=='juno')]), yend = sma30.coefs[3,1]+
                     sma30.coefs[3,2]*max(log10_beta[which(spec=='juno')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_beta[which(spec=='melpomene')]), xend = max(log10_beta[which(spec=='melpomene')]), 
                   y = sma30.coefs[4,1] + sma30.coefs[4,2]*min(log10_beta[which(spec=='melpomene')]), yend = sma30.coefs[4,1]+
                     sma30.coefs[4,2]*max(log10_beta[which(spec=='melpomene')])),colour="darkgrey",linewidth=1)+
  scale_colour_manual(values = c("#c22f20ca", "#1293b9ff", "#7dadb5ff","#cb8781ff"))
lobp30.plot.l

smatr.plot=ggarrange(lobp22.1.plot,lobp17.plot.l,lobp19.plot.l,lobp18.plot.l,lobp16.plot.l,lobp20.plot.l,
                     lobp8.plot.l,lobp8.1.plot,lobp9.plot.l,lobp7.plot.l,lobp6.plot.l,lobp10.plot.l,
                     lobp24.plot.l,lobp22.plot.l,lobp29.1.plot,lobp23.plot.l,lobp21.plot.l,lobp25.plot.l,
                     lobp13.plot.l,lobp12.plot.l,lobp14.plot.l,lobp15.1.plot,lobp11.plot.l, lobp15.plot.l,
                     lobp3.plot.l,lobp1.plot.l,lobp4.plot.l,lobp2.plot.l,lobp1.1.plot,lobp5.plot.l,
                     lobp29.plot.l,lobp27.plot.l,lobp30.plot.l,lobp28.plot.l,lobp26.plot.l,lobp36.1.plot,
                     ncol = 6, nrow = 6,common.legend = TRUE, legend="bottom")
smatr.plot


#repeat SMATR results for differences in clade (i.e. Heliconius vs outgroup)

sdata.clade=read.table(file="Dataset_LOB-FIN-STR_clade.txt", sep="\t", header=TRUE)
str(sdata.clade)

#filter for only young ones (removing species bias)
sdata.clade.young=filter(sdata.clade, age == "young")
str(sdata.clade.young)

Elev.test1.cl=sma(log10_gamma~log10_alpha_prime+clade, robust=T, data=sdata.clade.young)
summary(Elev.test1.cl)
Elev.test1.cl$gtr #to extract Wald stat.

Elev.test2.cl=sma(log10_gamma~log10_beta_prime+clade,robust=T,data=sdata.clade.young)
summary(Elev.test2.cl)
Elev.test2.cl$gtr

Elev.test3.cl=sma(log10_gamma~log10_alpha+clade,robust=T,data=sdata.clade.young)
summary(Elev.test3.cl)
Elev.test3.cl$gtr

Elev.test4.cl=sma(log10_gamma~log10_beta+clade,robust=T,data=sdata.clade.young)
summary(Elev.test4.cl)
Elev.test4.cl$gtr

Elev.test5.cl=sma(log10_gamma~log10_Y.lobe+clade,robust=T,data=sdata.clade.young)
summary(Elev.test5.cl)
Elev.test5.cl$gtr

Elev.test6.cl=sma(log10_alpha_prime~log10_gamma+clade,robust=T,data=sdata.clade.young)
summary(Elev.test6.cl)
Elev.test6.cl$gtr

Elev.test7.cl=sma(log10_alpha_prime~log10_beta_prime+clade,robust=T,data=sdata.clade.young)
summary(Elev.test7.cl)
Elev.test7.cl$gtr

Elev.test8.cl=sma(log10_alpha_prime~log10_alpha+clade,robust=T,data=sdata.clade.young)
summary(Elev.test8.cl)
Elev.test8.cl$gtr

Elev.test9.cl=sma(log10_alpha_prime~log10_beta+clade,robust=T,data=sdata.clade.young)
summary(Elev.test9.cl)
Elev.test9.cl$gtr

Elev.test10.cl=sma(log10_alpha_prime~log10_Y.lobe+clade,robust=T,data=sdata.clade.young)
summary(Elev.test10.cl)
Elev.test10.cl$gtr

Elev.test11.cl=sma(log10_beta_prime~log10_gamma+clade,robust=T,data=sdata.clade.young)
summary(Elev.test11.cl)
Elev.test11.cl$gtr

Elev.test12.cl=sma(log10_beta_prime~log10_alpha_prime+clade,robust=T,data=sdata.clade.young)
summary(Elev.test12.cl)
Elev.test12.cl$gtr

Elev.test13.cl=sma(log10_beta_prime~log10_alpha+clade,robust=T,data=sdata.clade.young)
summary(Elev.test13.cl)
Elev.test13.cl$gtr

Elev.test14.cl=sma(log10_beta_prime~log10_beta+clade,robust=T,data=sdata.clade.young)
summary(Elev.test14.cl)
Elev.test14.cl$gtr

Elev.test15.cl=sma(log10_beta_prime~log10_Y.lobe+clade,robust=T,data=sdata.clade.young)
summary(Elev.test15.cl)
Elev.test15.cl$gtr

Elev.test16.cl=sma(log10_alpha~log10_gamma+clade,robust=T,data=sdata.clade.young)
summary(Elev.test16.cl)
Elev.test16.cl$gtr

Elev.test17.cl=sma(log10_alpha~log10_alpha_prime+clade,robust=T,data=sdata.clade.young)
summary(Elev.test17.cl)
Elev.test17.cl$gtr

Elev.test18.cl=sma(log10_alpha~log10_beta_prime+clade,robust=T,data=sdata.clade.young)
summary(Elev.test18.cl)
Elev.test18.cl$gtr

Elev.test19.cl=sma(log10_alpha~log10_beta+clade,robust=T,data=sdata.clade.young)
summary(Elev.test19.cl)
Elev.test19.cl$gtr

Elev.test20.cl=sma(log10_alpha~log10_Y.lobe+clade,robust=T,data=sdata.clade.young)
summary(Elev.test20.cl)
Elev.test20.cl$gtr

Elev.test21.cl=sma(log10_beta~log10_gamma+clade,robust=T,data=sdata.clade.young)
summary(Elev.test21.cl)
Elev.test21.cl$gtr

Elev.test22.cl=sma(log10_beta~log10_alpha_prime+clade,robust=T,data=sdata.clade.young)
summary(Elev.test22.cl)
Elev.test22.cl$gtr

Elev.test23.cl=sma(log10_beta~log10_beta_prime+clade,robust=T,data=sdata.clade.young)
summary(Elev.test23.cl)
Elev.test23.cl$gtr

Elev.test24.cl=sma(log10_beta~log10_alpha+clade,robust=T,data=sdata.clade.young)
summary(Elev.test24.cl)
Elev.test24.cl$gtr

Elev.test25.cl=sma(log10_beta~log10_Y.lobe+clade,robust=T,data=sdata.clade.young)
summary(Elev.test25.cl)
Elev.test25.cl$gtr

Elev.test26.cl=sma(log10_Y.lobe~log10_gamma+clade,robust=T,data=sdata.clade.young)
summary(Elev.test26.cl)
Elev.test26.cl$gtr

Elev.test27.cl=sma(log10_Y.lobe~log10_alpha_prime+clade,robust=T,data=sdata.clade.young)
summary(Elev.test27.cl)
Elev.test27.cl$gtr

Elev.test28.cl=sma(log10_Y.lobe~log10_beta_prime+clade,robust=T,data=sdata.clade.young)
summary(Elev.test28.cl)
Elev.test28.cl$gtr

Elev.test29.cl=sma(log10_Y.lobe~log10_alpha+clade,robust=T,data=sdata.clade.young)
summary(Elev.test29.cl)
Elev.test29.cl$gtr

Elev.test30.cl=sma(log10_Y.lobe~log10_beta+clade,robust=T,data=sdata.clade.young)
summary(Elev.test30.cl)
Elev.test30.cl$gtr

####-----------------------####
#SMATr analysis of expanded structures vs non expanded, as posthoc to single pair-wise SMATrs.
Expansion_data=read.table(file="ExpandedParts_Data_.txt", sep="\t", header=TRUE)
str(Expansion_data)

Elev.testExp=sma(log10_Expanded~log10_NonExpanded+Species,multcomp=T, robust=T, data=Expansion_data)
summary(Elev.testExp)
Elev.testExp.grpsum=round(Elev.testExp$groupsummary[4], digits=3)
Elev.testExp.multcomp=round(Elev.testExp$multcompresult[ , c(3,4)], digits=3)
Elev.testExp$gtr

sma.coefs.exp=data.frame(coef(Elev.testExp))
sma.coefs.exp

#smatr analysis on clade level
Expansion_data.clade=read.table(file="ExpandedParts_Data_clade.txt", sep="\t", header=TRUE)
str(Expansion_data.clade)
Elev.testExp.clade=sma(log10_Expanded~log10_NonExpanded+clade, robust=T, data=Expansion_data.clade)
summary(Elev.testExp.clade)
Elev.testExp.clade$gtr

#plot
lobp.plot.exp=ggplot(Expansion_data,aes(x=log10_NonExpanded, y=log10_Expanded))+geom_point(aes(col=Species, size=0.1, alpha=0.5))+
  xlab(expression(bold(log[10]("not expanded"))))+ylab(expression(bold(log[10]("Expanded"))))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.border = element_blank(), panel.background = element_blank(),axis.line = element_line(color = 'black',size=0.8), legend.position = "none")
lobp.plot.exp
lobp.plot.exp2=lobp.plot.exp+
  geom_segment(aes(x = min(log10_NonExpanded[which(Species=='erato')]), xend = max(log10_NonExpanded[which(Species=='erato')]), 
                   y = sma.coefs.exp[1,1] + sma.coefs.exp[1,2]*min(log10_NonExpanded[which(Species=='erato')]), yend = sma.coefs.exp[1,1]+
                     sma.coefs.exp[1,2]*max(log10_NonExpanded[which(Species=='erato')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_NonExpanded[which(Species=='iulia')]), xend = max(log10_NonExpanded[which(Species=='iulia')]), 
                   y = sma.coefs.exp[2,1] + sma.coefs.exp[2,2]*min(log10_NonExpanded[which(Species=='iulia')]), yend = sma.coefs.exp[2,1]+
                     sma.coefs.exp[2,2]*max(log10_NonExpanded[which(Species=='iulia')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_NonExpanded[which(Species=='juno')]), xend = max(log10_NonExpanded[which(Species=='juno')]), 
                   y = sma.coefs.exp[3,1] + sma.coefs.exp[3,2]*min(log10_NonExpanded[which(Species=='juno')]), yend = sma.coefs.exp[3,1]+
                     sma.coefs.exp[3,2]*max(log10_NonExpanded[which(Species=='juno')])),colour="darkgrey",linewidth=1)+
  
  geom_segment(aes(x = min(log10_NonExpanded[which(Species=='melpomene')]), xend = max(log10_NonExpanded[which(Species=='melpomene')]), 
                   y = sma.coefs.exp[4,1] + sma.coefs.exp[4,2]*min(log10_NonExpanded[which(Species=='melpomene')]), yend = sma.coefs.exp[4,1]+
                     sma.coefs.exp[4,2]*max(log10_NonExpanded[which(Species=='melpomene')])),colour="darkgrey",linewidth=1)
lobp.plot.exp2
lobp.plot.exp3=lobp.plot.exp2+scale_colour_manual(values = c("#c22f20ca", "#1293b9ff", "#7dadb5ff","#cb8781ff"))
lobp.plot.exp3

####-----------------------####
#figure arrangement
fig.2.part=ggarrange(Rsq_CI.plot2,pca1.plot,lobp.plot.exp3,ncol = 3, nrow = 1)
fig.2
ncol = 6, nrow = 6,common.legend = TRUE, legend="bottom"

####-----------------------####
#cell numbers in GAD

setwd("C:/R-Analysis")
count.data=read.table(file="Cell_counts.txt", sep="\t", header=TRUE)
str(count.data)

count.lm=lm(count~spec,data=count.data)
summary(count.lm)
#0.002

count.plot=ggplot(count.data, aes(x=spec,y=count,fill=spec))
count.plot
count.plot2=count.plot+theme_minimal()+geom_dotplot(binaxis='y', stackdir='center',stackratio=1.5, dotsize=1.2)
count.plot2
count.plot3=count.plot2+scale_fill_manual(values = c("#1293b9ff", "#cb8781ff"))
count.plot3

#cell numbers in TH

setwd("C:/R-Analysis")
count.data.TH=read.table(file="Cell_counts_TH.txt", sep="\t", header=TRUE)
str(count.data.TH)

count.lm=lm(count~spec,data=count.data.TH)
summary(count.lm)
#0.913
count.data.TH$spec=factor(count.data.TH$spec,levels=
                           c("iulia","erato"))
count.plot=ggplot(count.data.TH, aes(x=spec,y=count,fill=spec))
count.plot
count.plot2=count.plot+theme_minimal()+geom_dotplot(binaxis='y', stackdir='center',stackratio=1.5, dotsize=1.2)
count.plot2
count.plot3=count.plot2+scale_fill_manual(values = c("#1293b9ff", "#c22f20ca"))
count.plot3