R/MVNH_functions.R

In lvmuyang/MVNH: Quantify and compare multivariate normal hypervolumes

## Functions for MVNH package
#library(MASS)
# db1 for dataframe1 db2 for dataframe2
# two random matrixes generated by 2 mutivariate normal distributions
#db1 = mvrnorm(n=50,mu=c(1,1),Sigma=matrix(c(1,0,0,1),nrow=2))
#db2 = mvrnorm(n=50,mu=c(2,2),Sigma=matrix(c(2,0.6,0.6,0.5),nrow=2))
#colnames(db1)= c("temp","prec")

#setwd("~/Desktop/Niche_overlap/MVNH")
# components of the determinant
MVNH_det = function(data=NULL,cov=NULL,var.names=NULL,cov.matrix=F,log=F){
  # fetch names from the data or covariance matrix
  if(!cov.matrix){
    if(is.null(var.names)){
      if(is.matrix(data) & !is.null(colnames(data))) var.name = colnames(data)
      if(is.matrix(data) & is.null(colnames(data))) var.name = paste0("variable",1:ncol(data))
      if(is.data.frame(data)) var.name = names(data)
    } else {
      var.name = var.names}
  } else{
    if(is.null(var.names)){
      if(is.matrix(cov) & !is.null(colnames(cov))) var.name = colnames(cov)
      if(is.matrix(cov) & is.null(colnames(cov))) var.name = paste0("variable",1:ncol(cov))
      if(is.data.frame(cov)) var.name = names(cov)
    } else {
      var.name = var.names}
  }

  if(cov.matrix){
    COV = cov
  } else {
    COV = cov(data)
  }
  s = diag(COV)
  rho = det(COV)/prod(s)
  det = c(det(COV),s,rho)
  names(det) = c("total",var.name,"cor")
  if(log) {
    return(log(det))}
  else{
    return(det)
  }
}
#MVNH_det(db1)

# components of the Bhattacharyya distance
MVNH_dissimilarity = function(db1=NULL,db2=NULL,mean1=NULL,mean2=NULL,cov1=NULL,cov2=NULL,metric="Bhattacharyya_distance",var.names=NULL,cov.matrix=F){
  if(!cov.matrix){
    S1 = cov(db1)
    S2 = cov(db2)
    S = (S1 + S2)/2
    u1 = apply(db1,2,mean)
    u2 = apply(db2,2,mean)
  } else{
    S1 = cov1
    S2 = cov2
    S = (S1 + S2)/2
    u1 = mean1
    u2 = mean2
  }
  # fetch names from the data or covariance matrix
  if(!cov.matrix){
    if(is.null(var.names)){
      if(is.matrix(db1) & !is.null(colnames(db1))) var.name = colnames(db1)
      if(is.matrix(db1) & is.null(colnames(db1))) var.name = paste0("variable",1:ncol(db1))
      if(is.data.frame(db1)) var.name = names(db1)
    } else {
      var.name = var.names}
  } else{
    if(is.null(var.names)){
      if(is.matrix(cov1) & !is.null(colnames(cov1))) var.name = colnames(cov1)
      if(is.matrix(cov1) & is.null(colnames(cov1))) var.name = paste0("variable",1:ncol(cov1))
      if(is.data.frame(cov1)) var.name = names(cov1)
    } else {
      var.name = var.names}
  }

  if(metric=="Bhattacharyya_distance"){
    # calculate the mahalnobis component
    bd_mah = as.vector((1/8)*t(u1-u2)%*%solve(S)%*%(u1-u2))
    bd_mah_1d = (1/8)*(u1-u2)^2/diag(S)
    bd_mah_cor = bd_mah - sum(bd_mah_1d)

    # calculate the determinant ratio
    bd_det = (1/2)*log(det(S)/(det(S1)*det(S2))^0.5)
    bd_det_1d = (1/2)*log(diag(S)/(diag(S1)*diag(S2))^0.5)
    bd_det_cor = bd_det-sum(bd_det_1d)

    # caculate the contribution of each variable and the correlation component
    bd = bd_mah+bd_det
    bd_1d = bd_mah_1d+bd_det_1d
    bd_cor = bd_mah_cor+bd_det_cor

    bd_mah_all = c(bd_mah,bd_mah_1d,bd_mah_cor)
    names(bd_mah_all) = c("total",var.name,"cor")
    bd_det_all = c(bd_det,bd_det_1d,bd_det_cor)
    names(bd_det_all) = c("total",var.name,"cor")
    bd_all = c(bd,bd_1d,bd_cor)
    names(bd_all) = c("total",var.name,"cor")

    return(list(Bhattacharyya_distance=bd_all,
                Mahalanobis_distance=bd_mah_all,
                Determinant_ratio=bd_det_all))
  }
  if(metric=="Pianka"){
    # calculate the mahalnobis component
    bd_mah = as.vector(t(u1-u2)%*%solve(S)%*%(u1-u2))
    bd_mah_1d = (u1-u2)^2/diag(S)
    bd_mah_cor = bd_mah - sum(bd_mah_1d)

    # calculate the determinant ratio
    bd_det = (1/2)*log(det(S)/(det(S1)*det(S2))^0.5)
    bd_det_1d = (1/2)*log(diag(S)/(diag(S1)*diag(S2))^0.5)
    bd_det_cor = bd_det-sum(bd_det_1d)

    # caculate the contribution of each variable and the correlation component
    bd = bd_mah+bd_det
    bd_1d = bd_mah_1d+bd_det_1d
    bd_cor = bd_mah_cor+bd_det_cor

    bd_mah_all = c(bd_mah,bd_mah_1d,bd_mah_cor)
    names(bd_mah_all) = c("total",var.name,"cor")
    bd_det_all = c(bd_det,bd_det_1d,bd_det_cor)
    names(bd_det_all) = c("total",var.name,"cor")
    bd_all = c(bd,bd_1d,bd_cor)
    names(bd_all) = c("total",var.name,"cor")

    return(list(Pianka=bd_all,
                Mahalanobis_distance=bd_mah_all,
                Determinant_ratio=bd_det_all))
  }
  if(metric=="MacArthur-Levins"){
    # calculate the mahalnobis component
    bd_mah = as.vector(t(u1-u2)%*%solve(S)%*%(u1-u2))
    bd_mah_1d = (u1-u2)^2/diag(S)
    bd_mah_cor = bd_mah - sum(bd_mah_1d)

    # calculate the determinant ratio
    bd_det = log(det(S)/det(S1))
    bd_det_1d = log(diag(S)/diag(S1))
    bd_det_cor = bd_det-sum(bd_det_1d)

    # caculate the contribution of each variable and the correlation component
    bd = bd_mah+bd_det
    bd_1d = bd_mah_1d+bd_det_1d
    bd_cor = bd_mah_cor+bd_det_cor

    bd_mah_all = c(bd_mah,bd_mah_1d,bd_mah_cor)
    names(bd_mah_all) = c("total",var.name,"cor")
    bd_det_all = c(bd_det,bd_det_1d,bd_det_cor)
    names(bd_det_all) = c("total",var.name,"cor")
    bd_all = c(bd,bd_1d,bd_cor)
    names(bd_all) = c("total",var.name,"cor")

    return(list(MacArthur_Levins=bd_all,
                Mahalanobis_distance=bd_mah_all,
                Determinant_ratio=bd_det_all))
  }
  if(metric=="Morisita"){
    # calculate the mahalnobis component
    bd_mah = as.vector(t(u1-u2)%*%solve(S)%*%(u1-u2))
    bd_mah_1d = (u1-u2)^2/diag(S)
    bd_mah_cor = bd_mah - sum(bd_mah_1d)

    # calculate the determinant ratio
    bd_det = log(det(S)^0.5*(det(S1)^0.5+det(S2)^0.5)/(2*(det(S1)*det(S2))^0.5))
    bd_det_1d = NA
    bd_det_cor = NA

    # caculate the contribution of each variable and the correlation component
    bd = bd_mah+bd_det
    bd_1d = NA
    bd_cor = NA

    bd_mah_all = c(bd_mah,bd_mah_1d,bd_mah_cor)
    names(bd_mah_all) = c("total",var.name,"cor")
    bd_det_all = c(bd_det)
    #names(bd_det_all) = c("total",var.name,"cor")
    bd_all = c(bd)
    #names(bd_all) = c("total",var.name,"cor")

    return(list(Morisita=bd_all,
                Mahalanobis_distance=bd_mah_all,
                Determinant_ratio=bd_det_all))
  }
}
#MVNH_dissimilarity(db1,db2,var.names = c("temp","prec"))
#I love Susie!