PGLMM.data.R

######################################
## Function that organizes the data so that PGLMM can be fit
## Created by Matthew Helmus
##
## attempt at fixing the function for use with larger data sets
## diag is replaced with Matrix:::Diagonal
######################################

require(Matrix)

fixed.pglmm.data <- function (modelflag = 1, sim.dat = NULL, samp = NULL, tree = NULL, 
          traits = NULL, env = NULL, Vcomp = NULL) 
{
  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, corr = 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, ])
  }
  Y <- samp
  U <- matrix(env, nrow = length(env), ncol = 1)
  X <- traits
  nsites <- dim(Y)[1]
  nspp <- dim(Y)[2]
  YY <- t(Y)
  YY <- as.matrix(as.vector(as.matrix(YY)))
  if (modelflag == 1) {
    Vfullspp <- kronecker(Diagonal(nsites), V)
    Vfullsite <- kronecker(Diagonal(nsites), Matrix(1, nspp, 
                                                nspp))
    VV <- list(Vfullspp = Vfullspp, Vfullsite = Vfullsite)
    XX <- kronecker(matrix(1, nsites, 1), Diagonal(nspp))
    return(list(YY = YY, VV = VV, XX = XX))
  }
  if (modelflag == 2) {
    u <- scale(U)
    U <- kronecker(u, Matrix(1, nspp, 1))
    Vfullspp <- kronecker(matrix(1, nsites, nsites), Diagonal(nspp))
    VfullsppV <- kronecker(matrix(1, nsites, nsites), V)
    VfullUCU <- Diagonal(as.vector(U)) %*% Vfullspp %*% Diagonal(as.vector(U))
    VfullUCUV <- Diagonal(as.vector(U)) %*% VfullsppV %*% Diagonal(as.vector(U))
    Vfullsite <- kronecker(Diagonal(nsites), Matrix(1, nspp, 
                                                nspp))
    VV <- list(VfullUCU = VfullUCU, VfullUCUV = VfullUCUV, 
               Vfullsite = Vfullsite)
    XXspp <- kronecker(matrix(1, nsites, 1), Diagonal(nspp))
    XX <- cbind(U, XXspp)
    YY <- as.vector(t(Y))
    return(list(YY = YY, VV = VV, XX = XX))
  }
  if (modelflag == 3) {
    u <- scale(U)
    U <- kronecker(u, Matrix(1, nspp, 1))
    Vfullsite <- kronecker(Diagonal(nsites), Matrix(1, nspp, 
                                                nspp))
    if (is.null(Vcomp)) {
      compscale <- 1
      Vcomp <- solve(V, Diagonal(nspp))
      Vcomp <- Vcomp/max(Vcomp)
      Vcomp <- compscale * Vcomp
    }
    Vfullcomp <- kronecker(Diagonal(nsites), Vcomp)
    VV <- list(Vfullcomp = Vfullcomp, Vfullsite = Vfullsite)
    XXspp <- kronecker(matrix(1, nsites, 1), Diagonal(nspp))
    XX <- cbind(XXspp, ((U %*% matrix(1, 1, nspp)) * XXspp))
    YY <- as.vector(t(Y))
    return(list(YY = YY, VV = VV, XX = XX))
  }
  if (modelflag == 4) {
    Vfulltrait <- kronecker(Diagonal(nsites), traits %*% t(traits))
    traitscale4 <- 100
    Vfulltrait <- traitscale4 * Vfulltrait
    Vfullsite <- kronecker(Diagonal(nsites), Matrix(1, nspp, 
                                                nspp))
    VV <- list(Vfulltrait = Vfulltrait, Vfullsite = Vfullsite)
    XXspp <- kronecker(matrix(1, nsites, 1), Diagonal(nspp))
    XX <- XXspp
    YY <- as.vector(t(Y))
    return(list(YY = YY, VV = VV, XX = XX))
  }
  if (modelflag == 5) {
    Vfulltrait <- kronecker(Diagonal(nsites), traits %*% t(traits))
    traitscale5 <- 10
    Vfulltrait <- traitscale5 * Vfulltrait
    Vfullspp <- kronecker(Diagonal(nsites), V)
    Vfullsite <- kronecker(Diagonal(nsites), Matrix(1, nspp, 
                                                nspp))
    VV <- list(Vfulltrait = Vfulltrait, Vfullspp = Vfullspp, 
               Vfullsite = Vfullsite)
    XXspp <- kronecker(matrix(1, nsites, 1), Diagonal(nspp))
    XX <- XXspp
    YY <- as.vector(t(Y))
    return(list(YY = YY, VV = VV, XX = XX))
  }
}