Simulates models used to test the PGLMM models

PGLMM.sim.R

# library(plotrix)
# library(ape)


# PGLMM_model_simulations.m

# Simulates models used to test the PGLMM models
# 23 June 2010 MATLAB CODE
# 6 September 2011 R CODE TRANSLATION - HELMUS

#Function Definition
PGLMM.sim <- function(tree,nsites=30,modelflag=1,figs=TRUE,second.env=TRUE,compscale = 1)
{
  if (!require(ape)) {stop("The 'ape' package is required")}
  if (!require(plotrix)) {stop("The 'ape' package is required")}
  
  bspp2 <- NULL
  Vcomp <- NULL
  envirogradflag2 <- 0

  if (is(tree)[1] == "phylo")
  {
    if (is.null(tree$edge.length))
    {
      tree <- compute.brlen(tree, 1)
    }
    V <- vcv.phylo(tree, cor = TRUE)
  } else {
    V <- tree
  }
   nspp<-dim(V)[1]
  # parameter values for each model
  if(modelflag==1 | modelflag==2)
  {
    Xscale <- 1
    Mscale <- .5
    Vscale1 <- 1
    Vscale2 <- 1
    b0scale <- .5

    envirogradflag1 <- 1
    if(second.env){envirogradflag2 <- 1} #envirogradflag2 <- 1
    elimsitesflag <- 1
    repulseflag <- 0
  }

  if(modelflag==3)
  {
    Xscale <- 1
    Mscale <- 0
    Vscale1 <- 1
    Vscale2 <- 1
    compscale <- compscale
    b0scale <- 0

    # repulsion matrix
    Vcomp <- solve(V,diag(nspp))
    Vcomp <- Vcomp/max(Vcomp)
    Vcomp <- compscale*Vcomp
    iDcomp <- t(chol(Vcomp))
    colnames(Vcomp)<-rownames(Vcomp)

    envirogradflag1 <- 1
    if(second.env){envirogradflag2 <- 1} # envirogradflag2=0;
    elimsitesflag <- 0
    repulseflag <- 1
  }

  if(modelflag==4)
  {
    Xscale <- 1
    Mscale <- .5

    Vscale1 <- .05
    Vscale2 <- 1

    b0scale <- .5

    envirogradflag1 <- 1
    if(second.env){envirogradflag2 <- 1} # envirogradflag2=0;
    elimsitesflag <- 1
    repulseflag <- 0
  }

  if(modelflag==5)
  {
    Xscale <- 1
    Mscale <- .5

    Vscale1 <- 1
    Vscale2 <- 1

    b0scale <- .5

    envirogradflag1 <- 1
    if(second.env){envirogradflag2 <- 1} #envirogradflag2=1

    elimsitesflag <- 1
    repulseflag <- 0
  }

  Vforsim <- V
  iD <- t(chol(Vforsim))

  # set up environmental gradient among sites

  mx <- t(as.matrix((-(nsites)/2):(nsites/2)))
  # number of sites
  m <- length(mx)


  #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  #% set up independent variables
  #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  #% environmental gradient 1
  if(envirogradflag1 == 1)
  {
    e <- iD %*% rnorm(nspp)
    e <- Vscale1*(e-mean(e))/sd(e)
    bspp1 <- e

    X <- 1/(1+exp(-(b0scale*array(1,c(m,1)) %*% rnorm(nspp) + t(mx) %*% t(e))))
    X <- Xscale*X
    Xsmooth <- X

    X1 <- diag(1-Mscale*runif(m)) %*% X
  }

  # environmental gradient 2
  if(envirogradflag2 == 1)
  {
    e <- iD %*% rnorm(nspp)
    e <- Vscale2*(e-mean(e))/sd(e)
    bspp2 <- e

    mx2 <- as.matrix(mx[sample(m)])
    X <- 1/(1+exp(-(b0scale*array(1,c(m,1)) %*% rnorm(nspp) + mx2 %*% t(e))))
    X <- Xscale*X
    Xsmooth <- Xsmooth*X

    X2 <- diag(1-Mscale*runif(m)) %*% X
    X <- X1*X2
  } else {
    X <- X1
  }

  # phylogenetic repulsion
  if(repulseflag == 1)
  {
    bcomp <- NULL
    for(i in 1:m)
    {
        bcomp <- cbind(bcomp, iDcomp %*% rnorm(nspp))
    }
    bcomp0 <- 0
    Xcomp <- exp(bcomp0+bcomp)/(1+exp(bcomp0+bcomp))
    X <- X*t(Xcomp)
  }

  #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  #% simulate data
  #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  #% convert distribution to presence/absence
  Y <- matrix(0,ncol=nspp,nrow=m)
  Y[matrix(runif(nspp*m),ncol=nspp) < X] <- 1
  colnames(Y)<-colnames(X)

  # eliminate sites with 0 spp
  pick <- (rowSums(Y)>0)
  Y <- Y[pick,]
  mx <- mx[pick]
  m <- length(mx)

  # eliminate sites with 1 spp
  if(elimsitesflag == 1){
    pick <- (rowSums(Y)>1)
    Y <- Y[pick,]
    mx <- mx[pick]
    m <- length(mx)
  }

  if(figs)
  {
    if (!require(plotrix))
    {
      stop("The 'plotrix' package is required to plot images from this function")
    }

    if(.Platform$OS.type == "unix") quartz() else windows()
    par(mfrow=c(5,1),las=1,mar=c(2, 4, 2, 2) - 0.1)
    matplot(Xsmooth,type="l",ylab="probability",main="Xsmooth")
    matplot(X,type="l",ylab="probability",main="X")
    hist(colSums(Y),main="spp per site")
    hist(rowSums(Y),main="sites per spp")
    plot(x=mx,y=rowSums(Y),main="SR across gradient",type="l")

    if(.Platform$OS.type == "unix") quartz() else windows()
    color2D.matplot(1-X/max(X),xlab="species",ylab="sites",main="probabilities")

    if(.Platform$OS.type == "unix") quartz() else windows()
    color2D.matplot(1-Y,xlab="species",ylab="sites",main="presence-absence")
  }

  return(list(Vphylo=V,Vcomp=Vcomp,Y=Y,X=X,u=mx,bspp1=bspp1,bspp2=bspp2))
}