structured.dmvn: compute multivariate normal densities for samples with group structure (e.g., individuals within species)

structured.dmvn

INSTRUCTIONS (see file below for script)

This is a slight modification of R's multivariate probability density function, dmvnorm, from the mvnorm package. The difference is that structured.dmvn allows means to vary by group (e.g., species, population, etc.). In contrast, dmvnorm evaluates all observations relative to a single, grand mean. The varying mean option is necessary when calculating log-likelihoods of observations in structured samples. 

INPUTS
1.	x: matrix of observations (indivs * traits)
2.	mean: matrix of means (indivs * traits)
3.	sigma: covariance matrix (traits * traits)

VALUES (returns)
1.	Log-likelihoods or likelihoods, depending on whether log=TRUE.

NOTES
1.	Other than the accommodation of varying means, function is same as dmvnorm.
2.	While means vary over the rows of the observations, sigma is shared by all.
3.	x may be a matrix or vector.

STEP-THROUGH
1.	Cholesky decomposition of covariance matrix.
2.	Compute densities or their logs.

structured.dmvn.r

structured.dmvn <- function (x, mean, sigma, log = TRUE) 
{
  #number of traits
  p <- ncol(x)

  #three tests for nonconformability
  if (ncol(mean) != p) 
      stop("mean and sigma have non-conforming size")
  
  if (ncol(sigma) != p)
      stop("x and sigma have non-conforming size")
  
  if (!isSymmetric(sigma, tol = sqrt(.Machine$double.eps), 
                     check.attributes = FALSE)) 
           stop("sigma must be a symmetric matrix")
  
  #test for whether sigma can be decomposed
  dec <- tryCatch(chol(sigma), error = function(e) e)
  if (inherits(dec, "error")) {
    x.is.mu <- colSums(t(x) != mean) == 0
    logretval <- rep.int(-Inf, nrow(x))
    logretval[x.is.mu] <- Inf
  }
  
  #calculate log likelihoods
  else {
    #changed t(x)-mean to t(x-mean) to accommodate varying means
    tmp <- backsolve(dec, t(x - mean), transpose = TRUE) 
    rss <- colSums(tmp^2)
    logretval <- -sum(log(diag(dec))) - 
      0.5 * p * log(2 * pi) - 0.5 * rss
  }
  
  #log likelhood or not
  names(logretval) <- rownames(x)
  if (log) 
    logretval
  else exp(logretval)
}