##Load raw data set
rd<-read.csv("rawdata.csv") 

##Prepare total lifed load data 
tll<-data.frame(tll=rd$Total.lifted.load,rd[,c(6:24)])

##Scale all varitbales
tll<-scale(tll)

##Fit the models
m.tll<-lm(tll~.,data = tll)

## pdredge the models
library(parallel)
library(MuMIn)
cores<-10
clusterType <- if(length(find.package("snow", quiet = TRUE))) "SOCK" else "PSOCK"
clust <- try(makeCluster(getOption("cl.cores", cores), type = clusterType))
clusterExport(clust, c("tll"), envir=environment())
system.time( dd<-pdredge(m.tll,rank = "AIC",cluster=clust) )
stopCluster(clust)
ddd <-subset(dd, delta < 6)
xx1<-model.avg(ddd)

# Create new data according to previos process, including important variables and their quadratic term
tll2<- data.frame(rd$Total.lifted.load, rd$Bill.length,rd$Bill.length^2,
                  rd$Wing.area,rd$Wing.area^2,
                  rd$Tarsus.length,rd$Tarsus.length^2,
                  rd$Flight.muscle.mass,rd$Flight.muscle.mass^2)
##Fit the models with quadratic term and interactions between important variables
m.tll.inter<-lm(tll~.*.,data = tll2)

## pdredge the models
cores<-10
clusterType <- if(length(find.package("snow", quiet = TRUE))) "SOCK" else "PSOCK"
clust <- try(makeCluster(getOption("cl.cores", cores), type = clusterType))
clusterExport(clust, c("tll2"), envir=environment())
system.time( dd<-pdredge(m.tll.inter,rank = "AIC",cluster=clust) )
stopCluster(clust)
ddd <-subset(dd, delta < 6)
xx2<-model.avg(ddd)

#Constructing SEM with the 4 most important variables
model1 <- "
# measurement model
Wing.size = ~ Wing.area
Flight.muscle.size = ~ chestto
Stomach.size = ~ Stomach.mass
Bill.size = ~  Bill.length
# regressions
TotalLiftedLoad ~ Wing.size + Flight.muscle.size +Stomach.size + Bill.size
"
#Fitting SEM with the 4 most important variables
fit <- sem(model, data=phy2)
summary(fit, standardized=TRUE,fit.measure = TRUE)

#Plotting SEM
semPaths(fit, "Standardized", "Estimates",rotation = 2,posCol=c("black","red"),levels = c(1,5,8),curve = 1.7,
         pastel = T,sizeMan = 15,sizeInt = 17,sizeMan2 = 15,sizeLat = 15,
         edge.label.cex =2, intercepts = T,layout = "tree",fade=F)