##### Supplementary code file 3/5   SET UP MODEL SPECIFICATIONS

##### M. Cardillo "Clarifying the relationship between body size and extinction risk in amphibians by complete mapping of model space"
##### marcel.cardillo@anu.edu.au
##### Apologies for the messiness and inelegance of the code!



library(caper)

##################################################################################
###### select Species, body_size_mm, aquatic, Lar, rangesize, latitude, Red.List.Ordinal, DD.ord.allLC, DD.ord.allEN, DD.ord.prop, hf, aridity, mpdq from the amphibio data frame

rows <- c(5,26,36,39,40,41,42,44,45,46,52,53,54)

##################################################################################
###### create a caper comparative dataset for each continent

compdat.africa <- comparative.data(phy=phy1,data=subset(amphibio, continent=="Africa",select=rows),names.col='Species')     # 504 rows
compdat.australia <- comparative.data(phy=phy1,data=subset(amphibio, continent=="Australia",select=rows),names.col='Species')     # 145 rows
compdat.eurasia <- comparative.data(phy=phy1,data=subset(amphibio, continent=="Eurasia",select=rows),names.col='Species')     # 491 rows
compdat.nthamerica <- comparative.data(phy=phy1,data=subset(amphibio, continent=="North America",select=rows),names.col='Species')     # 120 rows
compdat.sthamerica <- comparative.data(phy=phy1,data=subset(amphibio, continent=="South America",select=rows),names.col='Species')     # 1407 rows

##################################################################################
###### Set up model specifications:
# generic model is:
# test <- gls(response ~ bodysize + covariates,correlation=corStruct(1,phy=phy),data=compdat$data, na.action=na.omit,method='ML')

vars <- names(compdat.africa$data[c(2:5,10:12)]) # names of all the variables (apart from body size)
covariates <- unlist(lapply(1:length(vars), function(x) combn(vars, x, FUN = paste, collapse = "+")))
DD.species.ord <- c("omit","all.LC","all.EN","proportional")
phylogeny <- "phy1"

specs <- expand.grid(bodymass.values="length",covariates=c("none",covariates),DD.species=DD.species.ord,phylogeny=phylogeny,analysis=analysis.ord.phylo,IUCN='ordinal')[c(6:1)] # 10240 models


##################################################################################
###### Set up empty dataframe for results:

Models.pgls <- data.frame('fmla'=rep(NA,nrow(specs)),'fmla.null'=rep(NA,nrow(specs)),'correlation'=rep(NA,nrow(specs)),'variables'=rep(NA,nrow(specs)))
Models.pgls$AIC1 <- NA
Models.pgls$AIC0 <- NA
Models.pgls$BIC1 <- NA
Models.pgls$BIC0 <- NA
Models.pgls$int <-  NA
Models.pgls$slope <- NA
Models.pgls$slopeSE <- NA
Models.pgls$t <- NA
Models.pgls$p <- NA
Models.pgls$loglik1 <- NA
Models.pgls$loglik0 <- NA
Models.pgls$LR <- NA   # likelihood ratio test (Mcfadden's pseudo-rsq)
Models.pgls$RsqLik <- NA  # Ives' Rsqlik
Models.pgls$runtime <- NA

for (i in 1:nrow(specs)){
       # response
      if(specs$IUCN[i]=="ordinal" && specs$DD.species[i]=="omit") response <- "Red.List.Ordinal"
      if(specs$IUCN[i]=="ordinal" && specs$DD.species[i]=="all.LC") response <- "DD.ord.allLC"
      if(specs$IUCN[i]=="ordinal" && specs$DD.species[i]=="all.EN") response <- "DD.ord.allEN"
      if(specs$IUCN[i]=="ordinal" && specs$DD.species[i]=="proportional") response <- "DD.ord.prop"
      # bodysize
      if(specs$bodymass.values[i]=="length") bodysize <- "Body_size_mm"
      # covariates
      if(specs$covariates[i]=="none"){covariates <- NULL}else{covariates <- as.character(specs$covariates[i])}
       # model formula
      if(!is.null(covariates)){
           fmla <- paste(response,"~","log(",bodysize,")+",covariates,sep="")
           predictors <- c(bodysize,strsplit(covariates,split="\\+")[[1]])}
      if(is.null(covariates)){
           fmla <- paste(response,"~","log(",bodysize,")",sep="")
           predictors <- bodysize}
      # null formula       
      if(length(strsplit(outfile$variables[i],split=',')[[1]])< 3){ 
        fmla.null <- paste(response,"~1",sep="")}else{
        fmla.null <- paste(response,"~",specs$covariates[i],sep="")
        }   
    
      Models.pgls$fmla[i] <- fmla
      Models.pgls$fmla.null[i] <- fmla.null
      Models.pgls$variables[i] <- paste(c(response,predictors),collapse=",")
      }

Models.pgls[,c(6:18)] <- NA

##################################################################################
####### Duplicate Models.pgls and add body size^2 to the formulas, then rbind the two dataframes together

Models.pgls.2 <- Models.pgls
Models.pgls.2$fmla <- gsub("log(Body_size_mm)","log(Body_size_mm)+bmsq",Models.pgls.2$fmla,fixed=TRUE)
Models.pgls.2$variables <- gsub("Body_size_mm","Body_size_mm,bmsq",Models.pgls.2$variables,fixed=TRUE)
Models.pgls.2 <-  rbind(Models.pgls,Models.pgls.2)