Skip to content

Commit 2d7ea39

Browse files
pdpollypdpolly
authored
Add files via upload
1 parent ffe9746 commit 2d7ea39

2 files changed

Lines changed: 338 additions & 0 deletions

File tree

‎Goswami and Polly Modularity 2.0.m‎

Lines changed: 338 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,338 @@
1+
(* ::Package:: *)
2+
3+
(* This function prints the version number for the installed verison of this package. *)
4+
5+
ModularityVersion[]:=Print["Modularity for Mathematica 2.0\n(c) P. David Polly and A. Goswami, 20 Dec 2010\nUpdated 7 April 2018."];
6+
ModularityVersion[];
7+
8+
9+
(* Modularity Package, written by P.David Polly and Anjali Goswami 2010.
10+
11+
This package performs a simple analysis of modularity for geometric morphometric landmark data.
12+
13+
Three functions for landmarks are included, two of which calculate correlation matrices for multidimensional
14+
data (congruence coefficient and canonical correlation respectively) and one which uses those correlation
15+
matrices to perform a simple cluster analysis of modularity. An additional function, Magwene[],performs a graphical
16+
modeling analysis as described by Magwene (2001).The algorithms in the Magwene[] function do not work for singular
17+
correlation matrices,such as those derived from Procrustes superimposed landmarks.For further information see
18+
Goswami and Polly (2010). Functions for performing Procrustes analysis, matrix correlation analysis, Mantel's test, and
19+
subsampling analysis are also included.
20+
21+
Goswami,A. and P.D.Polly. 2010. Methods for studying morphological integration,modularity and covariance
22+
evolution.Pp.213-243 in J.Alroy and G.Hunt,eds., Quantitative Methods in Paleobiology. Paleontological
23+
Society Short Course, October 30th, 2010. The Paleontological Society Papers,Volume 16.
24+
25+
Magwene,P. 2001. New tools for studying integration and modularity. Evolution, 55:1734-1745.
26+
*)
27+
28+
29+
30+
(* This functions does a Procrustes superimposition of landmark coordinates followed by orthogonal projection into tangent space. The algorithm is the one
31+
presented by Rohlf & Slice, 1990, Syst. Zool. 39: 40-59 with tangent space projection from Rohlf, 1999, J. Class. 16: 197-223.
32+
33+
Usage: Procrustes[data, n, k]
34+
where data are the landmark coodiantes to be aligned, n is the number of landmarks, and k is the number of dimensions of each landmark (2 or 3).
35+
Aligned coordinates are returned with each shape in a single row, with the columns being the x,y (,z) coordinates of the n landmarks.
36+
37+
Updated 21 February 2010 to fix problem with 3D rotations.
38+
Updated 7 April 2018 to be consistent with current versions of Mathematica.
39+
*)
40+
Procrustes[data_,nlandmarks_,ndims_]:=Module[{l,II,PP,SS,x,y,u,w,v,hh,ResidSS,NewResidSS},
41+
l=Partition[Partition[Flatten[data],ndims],nlandmarks];
42+
II=IdentityMatrix[nlandmarks];
43+
PP=Table[N[1/nlandmarks],{nlandmarks},{nlandmarks}];
44+
SS=Table[N[Sqrt[Tr[(II-PP).l[[x]].Transpose[l[[x]]].(II-PP)]]],{x,Length[l]}];
45+
l=Table[((II-PP).l[[x]])/SS[[x]],{x,Length[l]}];
46+
y=Mean[l];
47+
ResidSS=Plus@@Flatten[Table[(Flatten[l[[x]]]-Flatten[y])^2,{x,Length[l]}]];
48+
While[True,
49+
For[x=1,x<=Length[l],x++,
50+
{u,w,v}=SingularValueDecomposition[Transpose[l[[x]]].y];
51+
hh=u.(w*Inverse[Abs[w]]).Transpose[v];
52+
l[[x]]=l[[x]].hh;
53+
54+
];
55+
y=Mean[l];
56+
NewResidSS=Plus@@Flatten[Table[(Flatten[l[[x]]]-Flatten[y])^2,{x,Length[l]}]];
57+
If[Abs[NewResidSS-ResidSS]<0.0001,Break[]];
58+
ResidSS=NewResidSS;
59+
];
60+
l=Partition[PrincipalComponents[Partition[Flatten[l],ndims]],nlandmarks];
61+
l=Partition[Flatten[l],ndims* nlandmarks];
62+
l=l.(IdentityMatrix[ndims nlandmarks]-Mean[l].Mean[l]);
63+
Return[l]]
64+
65+
66+
(* The CongruenceCoefficient[] function calculates a correlation matrix among a set of Procrustes superimposed
67+
2D or 3D landmarks the sum of the covariances between the two multidimensional landmarks over all specimens
68+
in the dataset divided by their pooled variance. Superimposed landmarks should be submitted in a three-
69+
dimensional array where the first dimension is the number of specimens, the second dimension is the number
70+
of landmarks,and the third dimension is the number of coordinates in each landmark. This function is used in
71+
the Modularity[] function below. Note that data must be partitioned after using the Procrustes function above before being
72+
entered in CongruenceCoefficient[] or CanonicalCorrelation[]. To partition data into 3D coordinates, use
73+
Partition[Partition[Flatten[superimposed],ndims],nlandmarks] as detailed in the user guide.
74+
*)
75+
CongruenceCoefficient[Superimposed_] := Module[ {i,j,N,stdevs},
76+
N = Table[dotcvmentry[Superimposed,i,j], {i,1,Dimensions[Superimposed][[2]]}, {j,1,Dimensions[Superimposed][[2]]}];
77+
stdevs=Sqrt[Tr[N,List]];
78+
Return[N/Table[Table[stdevs[[i]]*stdevs[[j]],{j,Length[stdevs]}],{i,Length[stdevs]}]]
79+
];
80+
81+
dotcvmentry[Superimposed_, col1_, col2_] := Module[ {i,n,s1,s2,p},
82+
n = 0;
83+
If[Dimensions[Superimposed][[3]]==2,
84+
s1 =s2= {0,0},s1 =s2= {0,0,0}];
85+
86+
For[i = 1, i <= Dimensions[Superimposed][[1]], i++,
87+
If[ !StringQ[Superimposed[[i,col1]]] && !StringQ[Superimposed[[i,col2]]],
88+
n++;
89+
s1 = s1 + Superimposed[[i,col1]];
90+
s2 = s2 + Superimposed[[i,col2]];
91+
]
92+
];
93+
If[n<=1,
94+
Print["there are too many missing data to covary columns ", col1, " and ", col2];
95+
Abort[];
96+
];
97+
s1 = s1/n;
98+
s2 = s2/n;
99+
100+
p = 0;
101+
For[i = 1, i <= Dimensions[Superimposed][[1]], i++,
102+
If[ !StringQ[Superimposed[[i,col1]]] && !StringQ[Superimposed[[i,col2]]],
103+
p = p + (Superimposed[[i,col1]] - s1) . (Superimposed[[i,col2]] - s2);
104+
];
105+
];
106+
p/(n-1)
107+
];
108+
109+
110+
111+
(* The CanonicalCorrelation[] function calculates a correlation matrix among a set of Procrustes superimposed
112+
2D or 3D landmarks as the correlation between the major axis scores of the respective landmarks. Superimposed
113+
landmarks should be submitted in a three-dimensional array where the first dimension is the number of
114+
specimens, the second dimsension is the number of landmarks,and the third dimension is the number of coordinates
115+
in each landmark. This function is used in the Modularity[] function below.
116+
117+
Updated 7 April 2018 to be consistent with current versions of Mathematica.
118+
*)
119+
120+
CanonicalCorrelation[Superimposed_]:=Module[{pcs,i},
121+
pcs=Transpose[Table[Transpose[PrincipalComponents[#-Mean[Superimposed[[1;;,i]]]&/@Superimposed[[1;;,i]]]][[1]],{i,Dimensions[Superimposed][[2]]}]];
122+
Return[Correlation[pcs]];
123+
];
124+
125+
126+
(* The Magwene[] function performs a graphical model analysis as described by Magwene 2001. "New tools for studying
127+
integration and modularity", Evolution, 55: 1734-1745. Input consists of a symmetric correlation matrix, a single
128+
integer for the sample size, and a vector of vertex labels. There are two optional variables, "Criterion" and
129+
"Threshold",which specify which matrix to use for edge labels and the threshold value for drawing a graph edge. By
130+
default the criterion is the p-value for edge exclusion deviance and the threshold is to draw edges for those
131+
whose probability is 0.95 or higher. Other criteria are "EED" for the edge exclusion deviance value itself, "ES"
132+
for edge strength, and "PC" for the partial correlation matrix. The edge labels are given from the same matrix
133+
specified by "Criterion". The function returns a graphical model for integration and modularity,
134+
the partial R-square values for each variable with respect to all others, the matrix of partial correlations
135+
among variables, the edge exclusion deviance matrix, the matrix of p-values for edge exclusion deviance, and the edge
136+
strength matrix.
137+
*)
138+
139+
Magwene[P_,SampSize_,VertexLabels_,Criterion_:"EEDP",Threshold_:0.95]:=Module[{\[CapitalOmega],\[Rho],i,j,RSquare,EED, EEDP, ES,GraphEdges,EdgeData},
140+
\[CapitalOmega] = Inverse[P];
141+
RSquare = (Tr[\[CapitalOmega],List]-1)/Tr[\[CapitalOmega],List];
142+
\[Rho] = Table[Table[-\[CapitalOmega][[i,j]]/(Sqrt[\[CapitalOmega][[i,i]]*\[CapitalOmega][[j,j]]]),{j,Length[P]}],{i,Length[P]}];
143+
Do[\[Rho][[i,i]]=1,{i,Length[\[Rho]]}];
144+
EED = Table[Table[(-SampSize)*Log[1-(\[Rho][[i,j]]^2)],{i,Length[\[Rho]]}],{j,Length[\[Rho]]}];
145+
EEDP=Chop[Table[Table[CDF[ChiSquareDistribution[1],EED[[i,j]]],{i,Length[\[Rho]]}],{j,Length[\[Rho]]}]];
146+
ES=Table[Table[(-0.5)*Log[1-(\[Rho][[i,j]]^2)],{i,Length[\[Rho]]}],{j,Length[\[Rho]]}];
147+
EdgeData=Switch[Criterion,
148+
"EEDP",EEDP,
149+
"EED",EED,
150+
"ES",ES,
151+
"PC",\[Rho]
152+
];
153+
GraphEdges=Flatten[Table[Table[{VertexLabels[[i]]->VertexLabels[[j]],EdgeData[[i,j]]},{j,i+1,Length[VertexLabels]}],{i,Length[VertexLabels]-1}],1];
154+
Return[{{GraphPlot[Select[GraphEdges,#[[2]]>Threshold&],VertexLabeling->True]},{"RSquare -> ",RSquare},{"Partial correlation matrix -> ",\[Rho]},{"Edge Exclusion Deviance Matrix ->",EED},{"Edge Exclusion Deviance P-values ->",EEDP},{"Edge Strength Matrix -> ",ES}}];
155+
]
156+
157+
158+
159+
(* The Modularity[] function does a simple analysis of modularity on geometric morphometric landmark data. The function requires the landmarks
160+
(formatted in a matrix where each row contains the x, y (, z) coordinates of all landmarks for a single specimen), the number of landmarks, the
161+
number of landmark dimensions (2 or 3), and labels for each landmark. The program returns results based on both the congruence coefficient and the
162+
canonical correlation (see Goswami and Polly, 2010). For each of the two coefficients the eigenvalues of the correlation matrix are returned in a
163+
barchart along with their standard deviation (see Pavlicev, Cheverud and Wagner, 2009), as well as a Ward's linkage dendrogram showing modular clusters
164+
among the landmarks.
165+
166+
Updated 7 April 2018 to be consistent with current versions of Mathematica.
167+
168+
*)
169+
<<HierarchicalClustering`
170+
Modularity[Landmarks_,NLands_,NDims_,Labels_,Randomization_:"False"]:=Module[{superimposed,residuals,P1,P2,EVStdev1,Plot1,EVStdev2,Plot2,
171+
NinetyFivePercentCutoff1,NinetyFivePercentCutoff2,Tree1,Tree2,EV1,EV2,BootEV1,BootEV2,MaxMinBootEV1,MaxMinBootEV2,MeanBootEV1,MeanBootEV2,
172+
Rows,Columns,dims,NumMods1,NumMods2,RelEVStdev1,RelEVStdev2},
173+
superimposed=Procrustes[Landmarks,NLands,NDims];
174+
superimposed=Partition[Partition[Flatten[superimposed],NDims],NLands];
175+
{Rows,Columns,dims}=Dimensions[superimposed];
176+
P1=CongruenceCoefficient[superimposed];
177+
P2=CanonicalCorrelation[superimposed];
178+
EV1=Chop[Eigenvalues[P1]];
179+
EVStdev1=Sqrt[Plus@@((Eigenvalues[P1]-Mean[Eigenvalues[P1]])^2)/MatrixRank[P1]];
180+
RelEVStdev1=Sqrt[Plus@@((Eigenvalues[P1]-Mean[Eigenvalues[P1]])^2)/MatrixRank[P1]]/Sqrt[MatrixRank[P1]];
181+
EV2=Chop[Eigenvalues[P2]];
182+
EVStdev2=Sqrt[Plus@@((Eigenvalues[P2]-Mean[Eigenvalues[P2]])^2)/MatrixRank[P2]];
183+
RelEVStdev2=Sqrt[Plus@@((Eigenvalues[P2]-Mean[Eigenvalues[P2]])^2)/MatrixRank[P2]]/Sqrt[MatrixRank[P2]];
184+
185+
186+
If[Randomization=="RandomizationTest",
187+
188+
BootEV1=Table[Eigenvalues[CongruenceCoefficient[Partition[superimposed[[#[[1]],#[[2]]]]&/@Flatten[ Table[Table[{RandomInteger[{1,Rows}],RandomInteger[{1,Columns}]},{Columns}],{Rows}],1],NLands]]],{100}];
189+
MaxMinBootEV1=Chop[Transpose[{Max[#]&/@Transpose[BootEV1],Min[#]&/@Transpose[BootEV1]}]];
190+
MeanBootEV1=Chop[Mean[#]&/@Transpose[BootEV1]];
191+
NumMods1=Length[Select[EV1-(#[[1]]&/@MaxMinBootEV1),Positive]];
192+
Plot1=Graphics[{
193+
{Gray,Table[Rectangle[{x-.45,0},{x+.45,EV1[[x]]}],{x,Length[EV1]}]},
194+
{Black,AbsoluteThickness[2],Dashed,Table[Line[{{x,MaxMinBootEV1[[x,2]]},{x,MaxMinBootEV1[[x,1]]}}],{x,Length[MaxMinBootEV1]}]},
195+
{Black,AbsolutePointSize[5],Table[Point[{x,MeanBootEV1[[x]]}],{x,Length[MeanBootEV1]}]}
196+
},PlotLabel->"Eigenvalues"];
197+
BootEV2=Table[Eigenvalues[CanonicalCorrelation[Partition[superimposed[[#[[1]],#[[2]]]]&/@Flatten[ Table[Table[{RandomInteger[{1,Rows}],RandomInteger[{1,Columns}]},{Columns}],{Rows}],1],NLands]]],{100}];
198+
MaxMinBootEV2=Chop[Transpose[{Max[#]&/@Transpose[BootEV2],Min[#]&/@Transpose[BootEV2]}]];
199+
MeanBootEV2=Chop[Mean[#]&/@Transpose[BootEV2]];
200+
NumMods2=Length[Select[EV2-(#[[1]]&/@MaxMinBootEV2),Positive]];
201+
Plot2=Graphics[{
202+
{Gray,Table[Rectangle[{x-.45,0},{x+.45,EV2[[x]]}],{x,Length[EV2]}]},
203+
{Black,AbsoluteThickness[2],Dashed,Table[Line[{{x,MaxMinBootEV2[[x,2]]},{x,MaxMinBootEV2[[x,1]]}}],{x,Length[MaxMinBootEV2]}]},
204+
{Black,AbsolutePointSize[5],Table[Point[{x,MeanBootEV2[[x]]}],{x,Length[MeanBootEV2]}]}
205+
},PlotLabel->"Eigenvalues"];
206+
Tree1=DendrogramPlot[DirectAgglomerate[Chop[1-Abs[P1]],Labels,Linkage->"Ward"],LeafLabels->(#&),HighlightLevel->If[NumMods1==0,None,NumMods1]];
207+
Tree2=DendrogramPlot[DirectAgglomerate[Chop[1-Abs[P2]],Labels,Linkage->"Ward"],LeafLabels->(#&),HighlightLevel->If[NumMods2==0,None,NumMods2]];
208+
Return[Style[Grid[{
209+
{Style["Modularity Results",FontWeight->Bold],""},
210+
{Style["Congruence Coefficient",FontSlant->Italic],""},
211+
{"Eignvalue SD: "<>ToString[EVStdev1]},
212+
{"Relative Eigenvalue SD: "<>ToString[RelEVStdev1]},
213+
{"Expectation of Random Rel. Eigenvalue SD: "<>ToString[Sqrt[1/Length[Landmarks]//N]]},
214+
{"Estimated Number of Modules: "<>ToString[NumMods1]},
215+
{Plot1,Tree1},
216+
{Style["Canonical Correlation",FontSlant->Italic],""},
217+
{"Eignvalue SD: "<>ToString[EVStdev2]},
218+
{"Relative Eigenvalue SD: "<>ToString[RelEVStdev2]},
219+
{"Expectation of Random Rel. Eigenvalue SD: "<>ToString[Sqrt[1/Length[Landmarks]//N]]},
220+
{"Estimated Number of Modules: "<>ToString[NumMods2]},
221+
{Plot2,Tree2}
222+
},Alignment->Left],FontFamily->"Helvetica"]
223+
]
224+
,
225+
226+
Plot1=BarChart[EV1,PlotLabel->"Eigenvalues"]; Plot2=BarChart[EV2,PlotLabel->"Eigvenvalues"];
227+
Tree1=DendrogramPlot[DirectAgglomerate[Chop[1-Abs[P1]],Labels,Linkage->"Ward"],LeafLabels->(#&)];
228+
Tree2=DendrogramPlot[DirectAgglomerate[Chop[1-Abs[P2]],Labels,Linkage->"Ward"],LeafLabels->(#&)];
229+
Return[Style[Grid[{
230+
{Style["Modularity Results",FontWeight->Bold],""},
231+
{Style["Congruence Coefficient",FontSlant->Italic],""},
232+
{"Eigenvalue SD: "<>ToString[EVStdev1]},
233+
{"Relative Eigenvalue SD: "<>ToString[RelEVStdev1]},
234+
{"Expectation of Random Rel. Eigenvalue SD: "<>ToString[Sqrt[1/Length[Landmarks]//N]]},
235+
{Plot1,Tree1},
236+
{Style["Canonical Correlation",FontSlant->Italic],""},
237+
{"Eigenvalue SD: "<>ToString[EVStdev2]},
238+
{"Relative Eigenvalue SD: "<>ToString[RelEVStdev2]},
239+
{"Expectation of Random Rel. Eigenvalue SD: "<>ToString[Sqrt[1/Length[Landmarks]//N]]},
240+
{Plot2,Tree2}
241+
},Alignment->Left],FontFamily->"Helvetica"]
242+
]
243+
244+
];
245+
246+
247+
248+
]
249+
250+
251+
252+
(*
253+
* Usage: mantel[A,B,n], where A and B are square matrices of the
254+
* same size, and n is an integer. Returns the matrix correlation
255+
* between the two original matrices, the fraction of times in the
256+
* n Monte Carlo randomization runs that the matrix correlation is
257+
* higher than a random correlation, and the probability that the
258+
* correlation is greater than random.
259+
*)
260+
261+
randperm[n_] := Module[ {l,t,r},
262+
l = Table[i, {i,1,n}];
263+
For [i = n, i > 1, i--,
264+
r = Random[Integer, {1,i}];
265+
t = l[[i]];
266+
l[[i]] = l[[r]];
267+
l[[r]] = t;
268+
];
269+
Return[l];
270+
];
271+
272+
randmperm[A_] := Module[ {n,p},
273+
If [!MatrixQ[A, NumberQ],
274+
Print["fatal error in mantel: bad matrix"];
275+
Abort[];
276+
];
277+
If [Dimensions[A][[1]] != Dimensions[A][[2]],
278+
Print["fatal error in mantel: matrix not square"];
279+
Abort[];
280+
];
281+
282+
n = Dimensions[A][[1]];
283+
p = randperm[n];
284+
Return[Table[ A[[p[[i]],p[[j]]]], {i,1,n}, {j,1,n} ]];
285+
];
286+
287+
mantelcor[l1_,l2_] := (l1.l2 - (Plus @@ l1) * (Plus @@ l2) / Length[l1]) /
288+
(Sqrt[l1.l1 - (Plus @@ l1) * (Plus @@ l1) / Length[l1]] *
289+
Sqrt[l2.l2 - (Plus @@ l2) * (Plus @@ l2) / Length[l2]]);
290+
291+
mantelmcor[A_,B_] := mantelcor[Flatten[A],Flatten[B]];
292+
293+
Mantel[A_,B_,n_] := Module[ {i,t,tt,npass},
294+
t = mantelmcor[A,B];
295+
npass = 0.0;
296+
For [i = 1, i <= n, i++,
297+
tt = mantelmcor[A,randmperm[B]];
298+
If [t >= tt, npass++]
299+
];
300+
Return[{t, npass/n, 1-npass/n}];
301+
];
302+
303+
304+
305+
(*Subsample[A,m]- subsamples data matrix A down to n rows*)
306+
RP[m_]:=Module[{rpx,rpi,rpj,rpt},
307+
rpx=Table[rpi,{rpi,1,m}];
308+
For[rpi=1,rpi<=m,rpi=rpi+1,rpj=Random[Integer,{1,m}];
309+
rpt=rpx[[rpi]];
310+
rpx[[rpi]]=rpx[[rpj]];
311+
rpx[[rpj]]=rpt;];
312+
Return[rpx];
313+
];
314+
315+
Subsample[A_,m_]:=Module[{rarx,rari,rary},
316+
rarx=RP[Length[A]];
317+
rary=Table[rarx[[rari]],{rari,1,m}];
318+
Return[A[[rary]]]
319+
];
320+
321+
322+
323+
(*subsamples matrix A from n-1 to n-m*)
324+
(*computes correlation coefficient with constant matrix B, can be same or different matrix*)
325+
(*k=number of landmarks, d=number of dimensions*)
326+
(*m=max.number of rows to remove*)
327+
(*t=no.of trials at each subsampling step*)
328+
329+
Diag[P_]:=Module[{i,j},Flatten[Table[Table[P[[i,j]],{j,i+1,Length[P]}],{i,Length[P]-1}]]]
330+
331+
SubsampleMatrix[A_,B_,nland_,ndims_,m_,t_]:=Module[{i,j,l},
332+
l=Length[A];
333+
Return[Table[
334+
Correlation[Diag[CongruenceCoefficient[Partition[Partition[Flatten[
335+
Procrustes[Subsample[A,l-i],nland,ndims]],ndims],nland]]],Diag[CongruenceCoefficient[Partition[Partition[Flatten[
336+
Procrustes[B,nland,ndims]],ndims],nland]]]],
337+
{i,1,m},{j,1,t}]]
338+
];

‎Modularity for Mathematica Users Guide 2.0.pdf‎

722 KB
Binary file not shown.

0 commit comments

Comments
 (0)