Function prepares data so that PGLMM can be fit

PGLMM.data.R

######################################
## Function that organizes the data so that PGLMM can be fit
## Created by Matthew Helmus
######################################

PGLMM.data<-function(modelflag=1,sim.dat=NULL,samp=NULL,tree=NULL,traits=NULL,env=NULL,Vcomp=NULL)
{ 
  if (!require(ape))
    {
      stop("The 'ape' package is required")
    }
  if(!is.null(sim.dat))
  {
    tree<-sim.dat$Vphylo
    Vcomp<-sim.dat$Vcomp
    samp<-sim.dat$Y
    traits<-sim.dat$bspp1
    env<-sim.dat$u
  }
  
  is.empty<-function(x){length(x) == 0}
  if(is.empty(samp))
  {
    stop("sample matrix (Y) is empty")
  }
  
  if (is(tree)[1] == "phylo")
  {
    if (is.null(tree$edge.length))
    {
      tree <- compute.brlen(tree, 1)
    }
    tree <- prune.sample(samp, tree)
    samp <- samp[, tree$tip.label]
    V <- vcv.phylo(tree, cor = TRUE)
    species <- colnames(samp)
    preval <- colSums(samp)/sum(sum(samp))
    species <- species[preval > 0]
    V <- V[species, species]
    Vcomp <- Vcomp[species, species]
    samp <- samp[, colnames(V)]
    traits <- as.matrix(traits[species,])
    

  } else {
    V <- tree
    species <- colnames(samp)
    preval <- colSums(samp)/sum(sum(samp))
    species <- species[preval > 0]
    V <- V[species, species]
    Vcomp <- Vcomp[species, species]
    samp <- samp[, colnames(V)]
    traits <- as.matrix(traits[species,])
    }
    
  # X = species independent variables (traits)
  # U = site independent variables (environment)
  # Y = site (rows) by species (columns) presence/absence matrix, the dependent variable (binary 0,1)
  # V = phylogenetic covariance matrix
  Y<-samp
  U<-matrix(env,nrow=length(env),ncol=1)
  X<-traits
  
  nsites<-dim(Y)[1]
  nspp<-dim(Y)[2]

  # create dependent variable vector
  YY=t(Y)
  YY=as.matrix(as.vector(as.matrix(YY)))

  if(modelflag==1)
  {
    # set up covariance matrices
    Vfullspp<-kronecker(diag(nsites),V)   #should be nspp*nsites in dimension
    Vfullsite<-kronecker(diag(nsites),matrix(1,nspp,nspp))

    #VV<-abind(Vfullspp,Vfullsite,along=3)  #could potentially also use the abind function in the abind library but I will see if we can get this to work instead
    VV<-list(Vfullspp=Vfullspp,Vfullsite=Vfullsite)
      
    # create independent variables
    XX<-kronecker(matrix(1,nsites,1),diag(nspp))
    return(list(YY=YY,VV=VV,XX=XX))
  }
 
  if(modelflag==2)
  {
    # set up covariance matrices
    u<-scale(U)
    U<-kronecker(u,matrix(1,nspp,1))
    
    Vfullspp<-kronecker(matrix(1,nsites,nsites),diag(nspp)) 
    VfullsppV<-kronecker(matrix(1,nsites,nsites),V)

    VfullUCU<-diag(as.vector(U)) %*% Vfullspp %*% diag(as.vector(U))
    VfullUCUV<-diag(as.vector(U)) %*% VfullsppV %*% diag(as.vector(U))

    Vfullsite<-kronecker(diag(nsites),matrix(1,nspp,nspp))

	  VV<-list(VfullUCU=VfullUCU,VfullUCUV=VfullUCUV,Vfullsite=Vfullsite)
    
    # create independent variables
    XXspp<-kronecker(matrix(1,nsites,1),diag(nspp))
    XX<-cbind(U,XXspp)

    # create dependent variable vector
    YY<-as.vector(t(Y))
    
    return(list(YY=YY,VV=VV,XX=XX))
  }
  
  if(modelflag == 3)
  {
    # set up covariance matrices
    u<-scale(U)
    U<-kronecker(u,matrix(1,nspp,1))
        
    #Vfullspp<-kronecker(diag(nsites),V)

    #VfullUCU<-diag(as.vector(U)) %*% Vfullspp %*% diag(as.vector(U))

    Vfullsite<-kronecker(diag(nsites),matrix(1,nspp,nspp))
    if(is.null(Vcomp))
    {
      compscale<-1
      # repulsion matrix
      Vcomp <- solve(V,diag(nspp))
      Vcomp <- Vcomp/max(Vcomp)
      Vcomp <- compscale*Vcomp
    }
    Vfullcomp<-kronecker(diag(nsites),Vcomp)

	  VV<-list(Vfullcomp=Vfullcomp,Vfullsite=Vfullsite)

    # create independent variables
    XXspp<-kronecker(matrix(1,nsites,1),diag(nspp))
    XX<-cbind(XXspp, ((U %*% matrix(1,1,nspp))*XXspp))

    # create dependent variable vector
    YY<-as.vector(t(Y))
    
    return(list(YY=YY,VV=VV,XX=XX))
  }
 
  if(modelflag == 4)
  {
    # set up covariance matrix for traits
    Vfulltrait <- kronecker(diag(nsites),traits %*% t(traits))
    traitscale4 <- 100
    Vfulltrait <- traitscale4*Vfulltrait
    # set up covariance matrices
    Vfullsite<-kronecker(diag(nsites),matrix(1,nspp,nspp))
 	  VV<-list(Vfulltrait=Vfulltrait,Vfullsite=Vfullsite)

    # create independent variables
    XXspp<-kronecker(matrix(1,nsites,1),diag(nspp))
    XX<-XXspp

    # create dependent variable vector
    YY<-as.vector(t(Y))
    return(list(YY=YY,VV=VV,XX=XX))
  }
  
  if(modelflag == 5)
  {  
    # set up covariance matrices
    Vfulltrait <- kronecker(diag(nsites),traits %*% t(traits))
    traitscale5 <- 10
    Vfulltrait <- traitscale5*Vfulltrait

    Vfullspp<-kronecker(diag(nsites),V)   #should be nspp*nsites in dimension

    Vfullsite<-kronecker(diag(nsites),matrix(1,nspp,nspp))
    
    VV<-list(Vfulltrait=Vfulltrait,Vfullspp=Vfullspp,Vfullsite=Vfullsite)
    
    # create independent variables
    XXspp<-kronecker(matrix(1,nsites,1),diag(nspp))
    XX<-XXspp

    # create dependent variable vector
    YY<-as.vector(t(Y))
    return(list(YY=YY,VV=VV,XX=XX))
  }
}

Dear Scott,

Have you managed to use these functions?
I have created a fork in which I try to fix some of the issues in the function (actually the one published in picante) that lead to memory being taken up and leading to system crashes.

Thank you very much