Added new function computeSAIR

computeSAIR computes the posterior distribution of the shared and idiosyncratic responses, as described in Ovaskainen and Abrego (manuscript).

R/computeSAIR.R

@@ -0,0 +1,108 @@
+#' @title computeSAIR
+#'
+#' @description Computes the shared and idiosyncratic responses to measured and latent predictors
+#'
+#' @param hM a fitted \code{Hmsc} model object
+#'
+#' @return
+#' returns the posterior distribution of the parameters mu.X2,V.XX,mu.Omega2,V.OmegaOmega,mu.tot,V.tot,s
+#' described Ovaskainen and Abrego (nmanuscript):
+#' How to synthesize a joint species distribution model into four numbers: shared and idiosyncratic responses of the species to measured and latent predictors
+#' 
+#' @details
+#' The shared and idiosyncratic responses are computed only for models without traits.
+#'
+#' @examples
+#' # Simulate a small dataset, fit Hmsc model to it, compute SAIR, and show posterior means 
+#'
+#' nc = 2
+#' ns = 5
+#' ny = 10
+#' mu = rnorm(n = nc)
+#' X = matrix(rnorm(n=nc*ny),nrow=ny)
+#' X[,1] = 1
+#' eps = matrix(rnorm(nc*ns),nrow=nc)
+#' L = matrix(rep(X%*%mu,ns),nrow=ny)  + X%*%eps
+#' Y = pnorm(L)
+#' m = Hmsc(Y = Y, XData = data.frame(env = X[,2]), distr = "probit")
+#' m = sampleMcmc(m,samples=100,transient=50,verbose = 0)
+#' SI = computeSAIR(m)
+#' colMeans(SI)
+#' 
+#' @importFrom stats cov
+#' @export
+
+
+computeSAIR = function(hM){
+
+  if(hM$nt>1) return(NULL)
+
+  if(hM$ncRRR==0){
+    switch(class(hM$X)[1L],
+           matrix = {
+             CX = cov(hM$X)
+           }, list = {
+             CX = lapply(hM$X, cov)
+           }
+    )
+  }
+
+  postList = poolMcmcChains(hM$postList)
+
+  n = length(postList)
+  SI = matrix(NA,ncol=7,nrow=n)
+  colnames(SI) = c("mu.X2","V.XX","mu.Omega2","V.OmegaOmega","mu.tot","V.tot","s")
+  for(i in 1:n){
+    post = postList[[i]]
+
+    if(hM$ncRRR>0){
+      XB=hM$XRRR%*%t(post$wRRR)
+      switch(class(hM$X)[1L],
+             matrix = {
+               CX = cov(cbind(hM$X,XB))
+             }, list = {
+               XX = hM$X
+               for(j in 1:length(XX)){
+                 XX[[j]] = cbind(XX[[j]],XB)
+               }
+               CX = lapply(XX, cov)
+             }
+      )
+    }    
+
+    mu = post$Gamma
+    V = post$V
+    switch(class(hM$X)[1L], matrix = {
+      mu.X2 = t(mu)%*%CX%*%mu
+      V.XX = sum(V*CX)
+    }, list = {
+      mu.X2 = 0
+      V.XX = 0
+      for(j in 1:hM$ns){
+        mu.X2 =  mu.X2 + t(mu)%*%CX[[j]]%*%mu
+        V.XX = V.XX + sum(V*CX[[j]])
+      }
+      mu.X2 =   mu.X2/hM$ns
+      V.XX =   V.XX/hM$ns
+    }
+    )
+    mu.Omega2 =  0
+    V.OmegaOmega = 0
+    if(hM$nr>0){
+      for(h in 1:hM$nr){
+        la = post$Lambda[[h]]
+        Omega = t(la)%*%la
+        Omega.diag = mean(diag(Omega))
+        diag(Omega) = NA
+        Omega.offdiag = mean(Omega, na.rm = TRUE)
+        mu.Omega2 =  mu.Omega2 + Omega.offdiag
+        V.OmegaOmega = V.OmegaOmega + (Omega.diag - Omega.offdiag)
+      }
+    }
+    mu.tot = mu.X2 + mu.Omega2
+    V.tot = V.XX + V.OmegaOmega
+    s = mu.tot/(mu.tot+V.tot)
+    SI[i,] = c(mu.X2,V.XX,mu.Omega2,V.OmegaOmega,mu.tot,V.tot,s)
+  }
+  return(SI)
+}