Loop Analysis in Mathematica

Loop analysis is an approach for calculating how press perturbations to a community might affect each species’ equilibrium abundance.  In simplified terms, it calculates the net effect of all the direct and indirect interactions between species given that only the presence and sign (+/-) of each pairwise interaction is known.

The following Mathematica commands calculate the classical adjoint (adjA, i.e. the net effects, a.k.a. the adjugate matrix, or the transpose of the cofactor matrix), as well as Dambacher et al.’s (2002) absolute feedback (T), and weighted-predictions (W) matrices given a so-called qualitatively-specified Community matrix.  See OSU’s Loop Group website for more info.  (Mathematica seems to use computer resources a bit more efficiently than Maple — for which Dambacher et al. provided code — especially for larger networks.)

Source: Novak, Wootton, Doak, Emmerson, Estes & Tinker (2011). Predicting community responses to perturbations in the face of imperfect knowledge and network complexity. Ecology, 92, 836-846.

Download: LoopAnalysis.m

(*Specify a community matrix*)

A = {{-1, 1, 1, 1},{-1, -1, 1, 1},{-1, -1, -1, 1},{-1, -1, -1, -1}}

(*Or, import community matrix from file in working directory*)


(*Create functions to compute matrix minors and matrix permanent*)

SetAttributes[ZD, Listable]; ZD[x_, y_]:= If[y == 0, 1, x/y];
Minor[m_List?MatrixQ,{i_Integer, j_Integer}]:=Abs[Drop[Transpose[Drop[Transpose[A],{j}]], {i}]]
Permanent[m_List]:= With[{v = Array[x, Length[m]]},Coefficient[Times @@ (m.v), Times @@ v]]

(*Calculate adjoint, absolute feedback, and weighted-predictions matrices *)

n = Length[A];
adjA = Inverse[-A]*Det[-A];
T = Outer[Permanent[Minor[Abs[A], {##}]] &, Sequence @@ Range /@ Dimensions[A], 1];
W = N[ZD[Abs[adjA], T]];
Export["AbsFeed.csv", T,"CSV"];

To source these Mathematica commands from within R, place the Mathematica Package (LoopAnalysis.m) into a working directory.  In R, export your qualitatively-specified community matrix (A) into the same working directly…

fout1<- file(paste('A.csv'),"w");

…and run one of the following lines (depending on your platform)…

shell('math < LoopAnalysis.m',mustWork=T) # in Windows
system('Mathkernel < LoopAnalysis.m',wait=T) # on Mac
system('math < LoopAnalysis.m',wait=T) # on server
…then import the results back into R and continue with your analyses.