Learning Interactions from MIcrobial Time Series (LIMITS)

Written by Charles K. Fisher (2014).

Feel free to use and modify this code how you wish, but please cite “Identifying Keystone Species in the Human Gut Microbiome from Metagenomic Timeseries using Sparse Linear Regression” by Fisher and Mehta.

RestrictedLeastSquares=Compile[{{inbag,_Real,2},{outbag,_Real,2},{test,_Integer,1}},Module[{X1,y1,B1,X2,y2,error},

(* This subroutine performs the linear regression and computes the error on the test set *)

(* 
"inbag" is a matrix containing the training set. One row of the matrix looks like {x_1(t) - u_1, ..., x_N(t) - u_N, ln x_i(t+1) - ln x_i(t) }, where u_i is the median abundance of species i. ;
"outbag" is a matrix containing the test set. It is in the same form as inbag.;
"test" is a vector of integers specifying which covariates are included in the regression. I.e. if species 1,2,3 are included then test = {1,2,3}. ;
*)

(* perform linear regression on training set *)
X1=Table[inbag[[k]][[#]]&/@test,{k,1,Length[inbag]}];
y1=Flatten[Take[inbag,All,-1]];
B1=PseudoInverse[X1].y1;

(* calculate prediction error on test set *)
X2=Table[outbag[[k]][[#]]&/@test,{k,1,Length[outbag]}];
y2=Flatten[Take[outbag,All,-1]];
error=Mean[(y2-X2.B1)^2]/Variance[y2];

(* export the error, the regression coefficients, and test in an array *)
{Table[error,{Length[test]}],B1,test}

],CompilationTarget->"C"];

LIMITS=Function[{data,n,i,thresh},Module[{BEST,percent,final,B,history,tmp2,excluded,included,testset,error=0.0,res,trainset,test,data1,data2,X1,y1,B1,X2,y2,r,prev},

(* This is the main function *)

(*
"data" is a matrix containing the data. One row of the matrix looks like {x_1(t) - u_1, ..., x_N(t) - u_N, ln x_i(t+1) - ln x_i(t) }, where u_i is the median abundance of species i. ;
"n" is the number of species ;
"i" is the species that is the dependent variable of the regression ;
"thresh" is the prediction error threshold; 
*)

r=100;(* the number of iterations used for Bagging ; *)
res={}; (* array for storing results *)

Do[

(* initialize covariates to be included/excluded from regression model *)
excluded=Delete[Range[n],i];
included={i};

(* randomly partiion data into training and test sets *)
trainset=Take[RandomSample[Range[Length[data]]],Ceiling[Length[data]/2.0]];
testset=Complement[Range[Length[data]],Union[trainset]];

data1=data[[trainset]];
data2=data[[testset]];

test=included;
BEST=RestrictedLeastSquares[data1,data2,test];

(* loop that adds covariates as long as prediction error decreases by thresh *) 
Do[
(* create list of performances for all regression models including an additional covariate *)
tmp2=Table[test=Append[included,excluded[[l]]];RestrictedLeastSquares[data1,data2,test],{l,1,Length[excluded]}];
(* sort the list so that prev[[1]] has the lowest prediction error *)
prev=Sort[tmp2,#1[[1]][[1]]<#2[[1]][[1]]&][[1]];
(* if prev[[1]] improves over the current best model by thresh, include it, otherwise exit loop *)
If[100.0*(prev[[1]][[1]]-BEST[[1]][[1]])/BEST[[1]][[1]]< -thresh, BEST=prev,Break[]];
(* redefine included/excluded covariates to account for including new covariate *)
included=IntegerPart[#]&/@BEST[[-1]];
excluded=Complement[excluded,included];
,{n-1}];

(* store final regression coefficients in res *)
B=Table[0.0,{n}];
(B[[#[[1]]]]=#[[2]])&/@Transpose[{BEST[[-1]],BEST[[2]]}];
AppendTo[res,B];

,{r}];

(* Bagging step: output median of res *)
Median[res]

]];

Example

RandomTimeSeries=Function[{P,sigma,L},

(* simulates a timeseries of length L for interaction matrix P and stochasticity sigma *)

Module[{n,x,tc},
n=Length[P];
tc=Table[0.0,{i,1,L+1},{j,1,n}];
x=RandomReal[LogNormalDistribution[0.0,sigma],{L+1,n}];
tc[[1]]=x[[1]];
Do[tc[[t+1]]=x[[t+1]]*tc[[t]]*Exp[P.Flatten[{1,tc[[t]]}]],{t,1,L}];
tc
]];

RandomAsymmetricInteractionMatrix[n_,seed_,sigma_,L_]:=Module[{eq,diag,const,PARAMS,foo, MaxEig = 0.9,goo,P,P1,k,l,tmp,tc,tc0},

(* creates a random interaction matrix that results in a stable timeseries *)

(* if seed for random number generation is provided, use it *)
If[seed≠0,SeedRandom[seed]];

(* randomly genrate equilibrium abundances *)
eq=RandomReal[LogNormalDistribution[0.0,0.1],n];

(* randomly construct a diagonal interaction matrix (species only interact with themselves) that produces the correct equilibrium abundances *)
diag=RandomReal[{-(1.0+MaxEig)/#,(MaxEig-1.0)/#}]&/@eq;
P=DiagonalMatrix[diag];
PARAMS=Transpose[Flatten[{{-P.eq},Transpose[P]},1]];

(* simulate a timeseries for this interaction matrix *)
tc0=RandomTimeSeries[PARAMS,sigma,L];

Do[

(* choose a position to try to add in a new interaction *)
{k,l}=RandomChoice[Position[P,0.0]];

(* randomly choose an interaction strength, which is less in absolute value than the self-interaction *)
tmp=(2.0*RandomReal[BetaDistribution[1.0,1.0]]-1.0)*Abs[diag[[k]]];

(* update interation matrix *)
P1=P;
P1[[k]][[l]]=tmp;
PARAMS=Transpose[Flatten[{{-P1.eq},Transpose[P1]},1]];

(* simulate a timeseries with the new interaction matrix *)
tc=Quiet[RandomTimeSeries[PARAMS,sigma,L]];

(* test for numerical stability, store interaction matrix if stable *)
If[MatrixQ[tc,MachineNumberQ[#]&]==False,
Break[];
,
P=P1;
tc0=tc;
];

,{t,1,L}];

(* ouptut interaction matrix and timeseries *)
{Transpose[Flatten[{{-P.eq},Transpose[P]},1]],tc0}
]

n=20;
thresh=5.0;

{PARAMS,tc}=RandomAsymmetricInteractionMatrix[n,0,0.05,400];
P=Take[PARAMS,All,{2,-1}];
ntc={(#/Total[#])&/@tc}; (* normalize timeseries to get relative abundances *)

ntc={tc};

u=Median[ntc[[1]]];
DATA=Table[Table[Table[Flatten[{ntc[[p]][[i]]-u,Log[ntc[[p]][[i+1]][[j]]]-Log[ntc[[p]][[i]][[j]]]}],{i,1,Dimensions[ntc][[2]]-1}],{p,1,Dimensions[ntc][[1]]}],{j,1,Dimensions[ntc][[3]]}];
BAGGING=Table[LIMITS[Flatten[DATA[[i]],1],n,i,thresh],{i,1,n}];

(* plot of timeseries *)
ListPlot[Transpose[ntc[[1]]],Joined->True,ImageSize->1000,AspectRatio->0.25,PlotRange->All]

GraphicsRow[{MatrixPlot[P],MatrixPlot[BAGGING]},ImageSize->500]

(* true interaction matrix *)
MatrixForm[P]

(* inferred interaction matrix scaled by population size *)
MatrixForm[Mean[u]*BAGGING]

Export["/Users/charlesfisher/Dropbox/learning_interactions/test_parameters.dat",P];

Export["/Users/charlesfisher/Dropbox/learning_interactions/test_data.dat",ntc[[1]]];

ntc[[1]]//MatrixForm