# The following R function (pgls.lam) computes a PGLS regression and finds the maximum likelihood value for ? using R’s 
# built-in optimising function ‘optim’. It uses the function ‘gls’ in the package ‘nlme’ to perform generalised least 
# squares regression with a supplied variance-covariance matrix (and hence requires that ‘nlme’ has been previously loaded). 
# We used the package ‘ape’ (Paradis, E, J. Claude and K. Strimmer. 2004. APE: analyses of phylogenetics and evolution in R 
# language. Bioinformatics 20:289-290.) to read in the phylogenetic tree, compute branch lengths according to 
# Grafen (1989) using the function ‘compute.brlen’ and to output the variance-covariance matrix for use in pgls.lam using 
# the function ‘vcv.phylo’ with the option ‘cor=TRUE’. Outputs from running pgls.lam include the estimated coefficients, 
# their standard errors and upper and lower 95% confidence intervals, the number of observations, number of estimated 
# parameters, log-likelihood, AIC, AICc (small sample version of AIC), the optimised value of ? and an indicator of 
# convergence (see the ‘optim’ function). An example call to pgls.lam is: 
#   model1 <- pgls.lam(model=log.rm ~ log.mass, VCV=mammal.vcv)

# pgls.lam fits phylogenetic least squares model assuming a Brownian motion model of evolution
# and optimises the value lambda - see Freckleton et al. (2002) American Naturalist 160: 712-726

  calc.lam <- function(lam,model,tree.weights,tree.cor) {
    m1 <- gls(model=model,cor=lam*tree.cor,weights=tree.weights,method="ML")
    return(logLik(m1))
  }

  pgls.lam <- function(model,VCV) {
    tree.weights <- varIdent(value=diag(VCV), fixed=T)
    tree.cor <- corSymm(VCV[lower.tri(VCV)],fixed=T)
    lam <- 1
    lmax <- optim(lam,fn=calc.lam,control=list(fnscale=-1),method="L-BFGS-B",lower=c(0),upper=c(1),model=model,tree.cor=tree.cor,tree.weights=tree.weights)
    model.lmax <<- gls(model=model,cor=lmax[[1]]*tree.cor,weights=tree.weights,method="ML")

    n.obs <- as.numeric(model.lmax$dims[1])
    loglik <- logLik(model.lmax)[1]
    n.obs <- attr(logLik(model.lmax),"nobs")
    n.par <- attr(logLik(model.lmax),"df") + 1       # add one for estimating lambda
    se <- sqrt(diag(model.lmax$varBeta))
    ucl <- coef(model.lmax) + qt(0.975,n.obs-n.par)*se
    lcl <- coef(model.lmax) - qt(0.975,n.obs-n.par)*se
    aic    <- -2*loglik + 2*n.par
    aicc   <- -2*loglik + 2*n.par + (2*n.par*(n.par+1))/(n.obs-n.par-1)
    pgls.lam.out <<- list(coef(model.lmax),se,ucl,lcl,n.obs,n.par,loglik,aic,aicc,lmax$par,lmax$convergence)
    names(pgls.lam.out) <<- c("summary","se","ucl","lcl","n.obs","n.par","logLik","AIC","AICc","lambda","converge")
    pgls.lam.out
  }