d<-read.table("Seabirddat_growth.txt",header=T, sep = "\t", dec = ".")
str(d)
attach(d)

## Load required packages
library(nlme)
library(caper)
library(ggplot2)

install.packages("ape")
library(ape)

## Load tree data from nexus file
trees <- read.nexus("output.nex")

## Have to create a species vector by hand for some reason...not able to connect the species of
## the phylogenetic trees to the species column in my data frame

row.names(d) <- c("Diomedea_exulans", "Diomedea_epomophora", "Thalassarche_chrysostoma", "Thalassarche_melanophrys", "Phoebetria_palpebrata", "Macronectes_giganteus", "Fulmarus_glacialis",
                  "Procellaria_aequinoctialis", "Puffinus_tenuirostris", "Puffinus_puffinus", "Oceanites_oceanicus", "Oceanodroma_leucorhoa", "Phaethon_rubricauda", 
                  "Pelecanus_occidentalis", "Fregata_magnificens", "Morus_bassanus", 
                  "Sula_dactylatra", "Sula_leucogaster", "Sula_sula", "Phalacrocorax_carbo", "Phalacrocorax_aristotelis", "Larus_dominicanus", "Larus_occidentalis", "Larus_argentatus", "Rissa_tridactyla", "Sterna_bergii",
                  "Sterna_sandvicensis", "Sterna_dougallii", "Sterna_hirundo", "Sterna_paradisaea", "Sterna_fuscata", "Anous_stolidus", "Anous_minutus", "Cepphus_columba", "Cepphus_grylle", "Uria_lomvia", "Alca_torda")

## Prune phylogenetic trees since the 3 species without growth rate information available were removed
class(trees)<-"multiPhylo"
pruned.trees<-lapply(trees,drop.tip,tip=c("Thalassarche_bulleri","Bulweria_bulwerii","Procellaria_parkinsoni"))
class(pruned.trees)<-"multiPhylo"

# Log-transform both predictor variable foraging range and dependent variable average weight
d$log.Forrang<-log(Forrang)
d$log.Weight_av<-log(Weight_av)

# More packages..
library(car)
library(AICcmodavg)

# Now trying to run the candidate models for PGLS model with growthrate
Cand.models = list()
niter = 100
for (i in 1:niter) {
  Cand.models[[i]] = gls(Maxgrowthrate ~ log.Forrang*log(Weight_av), data = d, method= "ML", 
                         correlation = corPagel(value=1, pruned.trees[[i]]))
}

## I get either error false convergence (8) (for lambda value = 1 for example) or an error in "eigen(val):infinite or missing values in 'x'" 
## (for lambda value = 0 for example)
## I don't really understand the problem because the datapoints look fine to me, except for the 
## few outliers but when removing these the same errors occur. There also no lack of datapoints, no overdispersion 


plot(d$log.Forrang, d$Maxgrowthrate, xlab = "log(foraging range (km))", ylab = "max. gr. rate (% ad. mass/d)",
     pch=16,cex=0.9,cex.lab=1.0,cex.axis=0.9)