Re: Compare matrices using Chi square for Exploratory
- To: mathgroup at smc.vnet.net
- Subject: [mg114942] Re: Compare matrices using Chi square for Exploratory
- From: Stuart Nettleton <Stuart.Nettleton at uts.edu.au>
- Date: Thu, 23 Dec 2010 03:55:40 -0500 (EST)
- References: <201012200539.AAA22748@smc.vnet.net>
For purposes of the record, the following m-file includes hypothesis tests
for both Maximum Likelihood and Variable variance. The former seems to be
highly unreliable in the contest of least squares, particularly if
recourse to PseudoInverse is necessary. The latter appears to be very
reliable, consistent with the communality of variables and therefore quite
usable. A calling file with some data and the final mfile are provided
below with the usual caveat of user beware and that any thoughts would be
much appreciated.
Many thanks to those who may have thought about this issue and warmest
seasons greetings to all!
Stuart
Calling program with some data considered to be very challenging:
Remove["Global`*"];
(*Set working directory where zigzagefa.m can be found*)
SetDirectory["path to your working directory"];
(*DATA1: Thurstone Box Problems from Thurstone 1940 1947 pp 141-44*)
tborigdata = {{3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5,
5, 5},
{2, 2, 3, 3, 3, 2, 2, 3, 3, 3, 4, 4, 4, 2, 2, 3, 3, 4, 4, 4},
{1, 2, 1, 2, 3, 1, 2, 1, 2, 3, 1, 2, 3, 1, 2, 2, 3, 1, 2, 3}};
(*26 variable test*)
tbvarsfn[{x_, y_, z_}] :=
{x, y, z, x y, x z, y z, x^2 y, x y^2,
x^2 z, x z^2, y^2 z, y z^2, x/y, y/x, x/z, y z/x, y/z, z/y,
2 x + 2 y, 2 x + 2 z, 2 y + 2 z, (x^2 + y^2)^(1/2), (x^2 + z^2)^(
1/2), (y^2 + z^2)^(1/2), x y z, (x^2 + y^2 + z^2)^(1/2)}
tbvars = N[Map[tbvarsfn, Transpose[tborigdata]]];
<< zigzagefa.m
?efa
?factorreports
{Z, \[CapitalLambda], \[CapitalPsi]result, factorscores,
frobeniusnormsqd,
maxiterations, \[ScriptCapitalH]1result, \[ScriptCapitalH]2result} \
= efa[tbvars, 3, Method -> "Triangular"];
factorreports[\[CapitalLambda], \[CapitalPsi] -> \
\[CapitalPsi]result,(*fs->factorscores,*)\[ScriptCapitalH]1 -> \
\[ScriptCapitalH]1result, \[ScriptCapitalH]2 -> \
\[ScriptCapitalH]2result,
desc -> "Unrotated"(*,vars->VariableNames,varslength->10*)]
MFILE zigzagefa.m
BeginPackage["zigzagefa`"]
(*If version before 8.0 then the BeginPackage line needs to be replaced=
with:*)
(*BeginPackage["zigzagefa`",{"MultivariateStatistics`"}]*)
efa::usage = "{Z, \[CapitalLambda], \[CapitalPsi], factor scores,
Frobenius \!\(\*SuperscriptBox[\(Norm\), \(2\)]\), maximum iterations,=
\[ScriptCapitalH]1 for a maximum liklihood test, \[ScriptCapitalH]2 for=
tests of variables} = efa[data, number of factors,
Method\[Rule]{Triangular(Default), Varimax, None}] is based on Unkel &=
Trendafilov's GEFALS (2010) for n>=p+k & Zig Zag EFA (2010) for n<p+k.=
Data is automatically corrected by setting non-numeric values to the
average for the variable. The EFA algorithm has been adjusted to set the=
factor maximums in the columns of \[CapitalLambda] to be positive at each=
iteration."
factorreports::usage="factorreports[\[CapitalLambda],
\[CapitalPsi]\[Rule]\[CapitalPsi] result (Default NA), fs\[Rule]factor=
scores result (Default none), \[ScriptCapitalH]1\[Rule]maximum liklihood=
test result (Default none), \[ScriptCapitalH]2\[Rule]test of variables=
result (Default none), vars\[Rule]variable names (Default is generic),=
varlength\[Rule]print length of variables (Default is 8),
desc\[Rule]optional prefix description (Default none)] prints the suite of=
reports."
varimaxTP::usage="varimaxTP[\[CapitalLambda], NormalizeVars\[Rule]{True,=
False(Default)}] performs an optimal varimax rotation. It operates on the=
factor loading matrix \[CapitalLambda] and returns \[CapitalLambda]vm.=
VarimaxTP uses the standard algorithm Kaiser(Psychometrika 1958) modified=
by Trevor Park in April 2002 from an original file by J.O.Ramsay. See
http://www.stat.ufl.edu/~tpark/Research/varimaxTP.m";
fixmissing::usage="fixmissing[data] sets any non-numeric data to the=
variable average. It is automatically used with efa"
Begin["`Private`"]
efa[data_,factors_,OptionsPattern[{Method->"Triangular"}]]:=Module[{fixed=
data,Z,n,p,k,B,RB,idx,F,U,\[CapitalLambda],\[CapitalPsi],f,f0,maxiter,m,U0,=
S0,V0,Q,uppert,Fbot,y1,yd,y2,=C3=9B,uW,dW,vW,\[CapitalLambda]negcolnum,temp=
2,fs,\[CapitalGamma],S,\[CapitalSigma],fk\[CapitalPsi],dof,\[Chi]2stat,\[Ch=
i]2pblty,\[Chi]2signif,\[Chi]2decision,\[Chi]2note,\[ScriptCapitalH]1,\[Scr=
iptCapitalH]2},
(*ZIGZAG EFA*)
(*Unkel S. & N.T. Trendfilov 2010 Zigzag exploratory factor analysis with=
more variables than observations*)
(*Z is mean centred, normalised to unit & variance adjusted*)
fixeddata=fixmissing[data];
Z=Transpose[Map[Normalize,Transpose[Map[#-Mean[fixeddata]&,fixeddata]]]];
{n,p}=Dimensions[Z];
k=factors; (*number of factors chosen exogenously*)
(*Initializing common factor scores F and unique factor scores U*)
If[n<p+k,
(*This is the original route of the Zig Zag paper with more variables than=
observations*)
{B,RB}=QRDecomposition[RandomReal[1,{p+k,n}]-0.5],
(*Use this for n>=p+k*)
{B,RB}=QRDecomposition[RandomReal[1,{n,n}]-0.5];
B=Take[B,All,p+k];
RB=Take[RB,p+k,p+k]
];
If[Min[Diagonal[RB]]<0,
idx=Pick[Range[Length[RB]],Thread[Diagonal[RB]<0]];
B[[idx,All]]=-B[[idx,All]]
];
F=Take[B,All,k];
U=Take[B,All,-p];
(*Initializing loadings \[CapitalLambda] and uniquenesses Psi*)
\[CapitalLambda]=Which[OptionValue[Method]=="Triangular",LowerTriangu=
larize[Transpose[Z].F],
OptionValue[Method]=="Varimax",varimaxTP[Transpose[Z].F],
OptionValue[Method]=="None",Transpose[Z].F,
True,Print["Error in Method Specification"];Abort[]
];
\[CapitalPsi]=DiagonalMatrix[Diagonal[Transpose[U].Z]];
\[CapitalGamma]=F.Transpose[\[CapitalLambda]]+U.\[CapitalPsi];
f=Norm[Z-\[CapitalGamma],"Frobenius"]^2;(*% Objective function*)
f0=f-1;
(*Main loop for ALS algorithm*)
maxiter=0;
While[Abs[f-f0]>10^-6 ,(*% Convergence criterion*)
(*Solve Procrustes problem to Method F*)
If[n<p+k,
(*This is the original route of the Zig Zag paper with more variables than=
observations*)
m=(Z-U.\[CapitalPsi]).\[CapitalLambda];
{U0,S0,V0}=SingularValueDecomposition[m,Length[SingularValueList[m]]];
F=U0.Transpose[V0];
(*Solve reduced Procrustes-like problem for =C3=9B,such that*)
(*=C3=9B.Transpose[=C3=9B]=IdentityMatrix[n-k] and hence
F.Transpose[F]+U.Transpose[U]=IdentityMatrix[n]*)
{Q,uppert}=QRDecomposition[Map[PadRight[#,n]&,F]];
Q=Transpose[Q];
Fbot=Take[Q,All,{k+1,n}];
m=\[CapitalPsi].Transpose[Z].Fbot;
{y1,yd,y2}=SingularValueDecomposition[m,Length[SingularValueList[m]]];
=C3=9B=y2.Transpose[y1];
(*Method U*)
U=Fbot.=C3=9B,
(*From GEFALS for n>=p+k*)
m=Z.Join[\[CapitalLambda],\[CapitalPsi],2];
{uW,dW,vW}=SingularValueDecomposition[m,Length[SingularValueList[m]]];
B=uW.Transpose[vW];
F=Take[B,All,k];
U=Take[B,All,{k+1,k+p}]
];
(**)
(*Method \[CapitalLambda] and Psi*)
\[CapitalLambda]=Which[OptionValue[Method]=="Triangular",LowerTriangu=
larize[Transpose[Z].F],
OptionValue[Method]=="Varimax",varimaxTP[Transpose[Z].F],
OptionValue[Method]=="None",Transpose[Z].F,
True,Print["Error in Method Specification"];Abort[]
];
(*Set the maximum value in each column to be positive and correspondingly=
adjust F*)
\[CapitalLambda]negcolnum=Position[Map[Sign,Map[\[CapitalLambda][[Orderin=
g[Transpose[Abs[\[CapitalLambda]]][[#]],-1][[1]],#]]&,Range[k]]],-1];
If[Length[\[CapitalLambda]negcolnum]>0,
\[CapitalLambda]=Transpose[MapAt[Minus,Transpose[\[CapitalLambda]],\[Capi=
talLambda]negcolnum]];
F=Transpose[MapAt[Minus,Transpose[F],\[CapitalLambda]negcolnum]]
];
\[CapitalPsi]=DiagonalMatrix[Diagonal[Transpose[U].Z]];
(*Evaluating the objective function*)
f0=f;
\[CapitalGamma]=F.Transpose[\[CapitalLambda]]+U.\[CapitalPsi];
f=Norm[Z-\[CapitalGamma],"Frobenius"]^2;
If[Mod[maxiter,100]==0,NotebookDelete[temp2];temp2=PrintTemporary["It=
eration
"<>ToString[maxiter]]];
maxiter++
(*End of main loop*)
];
(*factor scores*)
fs=Z.\[CapitalLambda];
\[CapitalSigma]=\[CapitalLambda].Transpose[\[CapitalLambda]]+\[CapitalPsi=
].Transpose[\[CapitalPsi]];
S=Covariance[Z];
(*Maximum Likelihood Hypothesis*)
(*Lawley, D.N. & Maxwell, M.A. Factor Analysis as a Statistical Method,=
Butterworths, London, 1971, Section 4.4 Tests of Hypotheses pp34-8 &
92-93*)
(*Note S.Inverse[\[CapitalSigma]]==p & L&M p35 provides a more stable=
alternative for
S.Inverse[\[CapitalSigma]]=IdentityMatrix[p]-(\[CapitalSigma]-S).PseudoIn=
verse[\[CapitalPsi]^2]]*)
(*The following is the same for the L&M model & R model which uses
fk\[CapitalPsi]=Tr[S.Inverse[\[CapitalSigma]]]\[Minus]Log[Det[S.Inverse[\=
[CapitalSigma]]]]\[Minus]p*)
fk\[CapitalPsi]=Quiet[Check[
If[MatrixQ[Inverse[\[CapitalSigma]]],Re[Log[Det[\[CapitalSigma]]]-Log[Det[S=
]]+Tr[S.Inverse[\[CapitalSigma]]]-p]],
(*Re[Log[Det[\[CapitalSigma]]]-Log[Det[S]]+Tr[S.PseudoInverse[\[CapitalSigm=
a]]]-p]*)
(*Form 1 Using direct PseudoInverse*)
Re[Log[Det[\[CapitalSigma]]]-Log[Det[S]]+Tr[IdentityMatrix[p]-(\[CapitalSig=
ma]-S).PseudoInverse[\[CapitalPsi]^2]]-p]
(*Form 2 from L&M p35 Eq 4.27*)
(*Re[Tr[S.PseudoInverse[\[CapitalSigma]]]-Log[Det[S.PseudoInverse[\[Capital=
Sigma]]]]-p]*)
(*Form 3 from R package factanal, see Psych package documentation*)
(*Re[Log[Det[\[CapitalSigma]]]-Log[Det[S]]]*)(*Form 4 based on L&M and=
Grolush p .158 quoting L&M eqn 7.15 assuming that
S.Inverse[\[CapitalSigma]]==p*)
]];
dof=((p-k)^2-(p+k))/2;
\[Chi]2stat=(n\[Minus]1\[Minus](2 p+5)/6\[Minus]2 k/3) fk\[CapitalPsi] ;
\[Chi]2pblty={0.01,0.05,0.10};
\[Chi]2signif=Map[InverseCDF[ChiSquareDistribution[dof],#]&,\[Chi]2pblty];
\[Chi]2decision=Map[\[Chi]2stat>#&,\[Chi]2signif]/.{True->"Reject
\[ScriptCapitalH]1",False->"Accept \[ScriptCapitalH]1"};
\[Chi]2note=Quiet[Check[If[MatrixQ[Inverse[\[CapitalSigma]]],{"using=
direct Inverse of \[CapitalSigma] (Pseudo Inverse not required)"}],{"using=
Pseudo Inverse of \[CapitalSigma]"}]];
\[ScriptCapitalH]1={\[Chi]2stat,dof,\[Chi]2pblty,Transpose[{\[Chi]2signif=
,\[Chi]2decision}],\[Chi]2note};
\[ScriptCapitalH]2=Map[VarianceEquivalenceTest[{\[CapitalGamma][[All,#]],=
Z[[All,#]]},"HypothesisTestData"]&,Range[p]];
Print["Goodness of fit (Frobenius \!\(\*SuperscriptBox[\(Norm\),
\(2\)]\)): ",f];
Return[{Z,\[CapitalLambda],Diagonal[\[CapitalPsi]],fs,f,maxiter,\[ScriptCap=
italH]1,\[ScriptCapitalH]2}]
];
fixmissing[data_]:=Module[{means,missing},
means=Map[Mean[Select[Transpose[data][[#]],NumberQ]]&,Range[Dimensions[da=
ta][[2]]]];
missing=Position[Map[Map[NumericQ,#]&,data],False];
Return[ReplacePart[data,Thread[missing->means[[missing[[All,2]]]]]]]
];
tableout[data_,tablename_,rownames_,colnames_,OptionsPattern[{rounding->0.0=
01,units->""}]]:=Module[{upper,lower},
If[Length[Dimensions[data]]==1,
upper={{Text[tablename],SpanFromLeft},{Text[OptionValue[units]],colnames}=
};
lower=Transpose[{rownames,If[VectorQ[data,NumberQ],Round[data,OptionValue=
[rounding]],data]}],
upper={Join[{Text[tablename]},Table[SpanFromLeft,{j,Dimensions[data][[2]]=
}]],Join[{Text[OptionValue[units]]},colnames]};
lower=Transpose[Prepend[Map[If[VectorQ[#,NumberQ],Round[#,OptionValue[rou=
nding]],#]&,Transpose[data]],rownames]]
];
Print[Grid[Join[upper,lower],Background->{None,{Lighter[Yellow,.9],{White,L=
ighter[Blend[{Blue,Green}],.8]}}},
Dividers->{{Darker[Gray,.6],{Lighter[Gray,.5]},Darker[Gray,.6]},{Darker[Gra=
y,.6],Darker[Gray,.6],{False},Darker[Gray,.6]}},
Alignment->{{Left,Center,{Center}}},Frame->Darker[Gray,.6],ItemStyle->14,Sp=
acings->{Automatic,.8}]]
];
factorreports[\[CapitalLambda]_,OptionsPattern[{\[CapitalPsi]->{},fs->{},\[=
ScriptCapitalH]1->{},\[ScriptCapitalH]2->{},vars->{},varslength->8,desc->""=
}]]:=Module[{\[CapitalPsi]2,\[CapitalLambda]1,\[CapitalLambda]2,\[Capital=
Lambda]3,\[CapitalLambda]4,\[CapitalLambda]5,p,k,facnames,varnames,signames=
,rinames,obsnames},
{p,k}=Dimensions[\[CapitalLambda]];
facnames=Table["Fac"<>ToString[i],{i,k}];
varnames=Table["Var"<>ToString[i],{i,p}];
If[Length[OptionValue[vars]]>0,
If[Length[OptionValue[vars]]!=p,
Print["Warning, list of variable names different to the number of
variables, continuing with generic variable names ...."],
varnames=Map[StringTake[ToString[OptionValue[vars][[#]]],Min[StringLength=
[ToString[OptionValue[vars][[#]]]],OptionValue[varslength]]]&,Range[Length[=
OptionValue[vars]]]]
]];
signames=Table["Sig"<>ToString[i],{i,k}];
rinames=Table["RI"<>ToString[i],{i,k}];
obsnames=Table["Obs"<>ToString[i],{i,Dimensions[OptionValue[fs]][[1]]}];
If[Length[OptionValue[\[CapitalPsi]]]==0,\[CapitalPsi]2=Table["NA",{i=
,p}],\[CapitalPsi]2=OptionValue[\[CapitalPsi]]^2];
\[CapitalLambda]1=\[CapitalLambda]^2;
\[CapitalLambda]2=Total[\[CapitalLambda]1]/p;
\[CapitalLambda]3=Join[\[CapitalLambda]1,{\[CapitalLambda]2},{Total[\[Cap=
italLambda]1]/Total[\[CapitalLambda]1,2]}];
\[CapitalLambda]4=MapThread[Append,{\[CapitalLambda]3,Total[\[CapitalLamb=
da]3,{2}]}];
\[CapitalLambda]5=Join[\[CapitalLambda]4,Map[\[CapitalLambda]3[[#]]/Total=
[\[CapitalLambda]3,{2}][[#]]&,Range[Length[\[CapitalLambda]3]]],2];
Print[Column[{
tableout[Join[MapThread[Append,{\[CapitalLambda],\[CapitalPsi]2}],Map[Cases=
[#,a_->Which[Abs[a]>=0.9,">.9",
Abs[a]>=0.7,">.7",Abs[a]>=
0.6,">.6",Abs[a]>=0.4,">.4",Abs[a]>=0.25,">.25",True,""]]&,\[CapitalLam=
bda]],2],
OptionValue[desc]<>": Factor Loadings (Fac), Residuals
(\!\(\*SuperscriptBox[\(\[CapitalPsi]\), \(2\)]\)) &
Significance",varnames,
Flatten[{facnames,{"\[CapitalPsi]^2"},signames}],rounding->0.01,units->"R"],
If[Length[OptionValue[\[ScriptCapitalH]1]]==0,"",
If[OptionValue[\[ScriptCapitalH]1][[1]]<0,
" \!\(\*
StyleBox[\"NOTE\",\nFontWeight->\"Bold\"]\)\!\(\*
StyleBox[\" \",\nFontWeight->\"Bold\"]\)\!\(\*
StyleBox[\"\[ScriptCapitalH]1\",\nFontWeight->\"Bold\"]\)\!\(\*
StyleBox[\" \",\nFontWeight->\"Bold\"]\)\!\(\*
StyleBox[\"Maximum\",\nFontWeight->\"Bold\"]\)\!\(\*
StyleBox[\" \",\nFontWeight->\"Bold\"]\)\!\(\*
StyleBox[\"Likelihood\",\nFontWeight->\"Bold\"]\)\!\(\*
StyleBox[\" \",\nFontWeight->\"Bold\"]\)\!\(\*
StyleBox[\"\[Chi]2\",\nFontWeight->\"Bold\"]\)\!\(\*
StyleBox[\" \",\nFontWeight->\"Bold\"]\)\!\(\*
StyleBox[\"test\",\nFontWeight->\"Bold\"]\):
The reason that this note has been printed is that in this particular=
extraction,
\[Chi]2stat is negative (" <>
ToString[Round[OptionValue[\[ScriptCapitalH]1][[1]],1]] <> " " <>
OptionValue[\[ScriptCapitalH]1][[5]] <> "), which is not an
unsurprising outcome for this calculation. \[Chi]2stat may be positive if=
the factors
are extracted again. The \[ScriptCapitalH]1 Maximum Likelihood \[Chi]2=
test can be highly unreliable
for random start least squares, particularly if the Pseudo Inverse of=
\[CapitalSigma] is required.
This test is provided for interest only and should not be used as a=
basis to
accept or reject. It is recommeded that the \[ScriptCapitalH]2 variables=
test and Communalities
be relied upon instead.",
tableout[OptionValue[\[ScriptCapitalH]1][[4]],
OptionValue[desc]<> ": " <> "This \[ScriptCapitalH]1 Maximum Likelihood=
\[Chi]2 test can be highly unstable for random start least squares,
particularly if the Pseudo Inverse of \[CapitalSigma] is identified as=
being required (see \[Chi]2stat below).
This test is provided for interest only and should not be used as a basis=
to accept or reject.
It is recommeded that the \[ScriptCapitalH]2 variables test and
Communalities be relied upon used instead.
Test of \[ScriptCapitalH]1 that " <> ToString[k] <> " factors is
sufficient:
\[Chi]2stat: " <> ToString[Round[OptionValue[\[ScriptCapitalH]1][[1]],1]]=
<> " " <> OptionValue[\[ScriptCapitalH]1][[5]] <> ".
Degrees of freedom: " <> ToString[OptionValue[\[ScriptCapitalH]1][[2]]]=
<> ".
Reject \[ScriptCapitalH]1 if \[Chi]2statistic exceeds the significance=
level for \[Chi]2 corresponding to a significance probability.",
Map[ToString,OptionValue[\[ScriptCapitalH]1][[3]]],{"Significance
Level","Conclusion"},rounding->1,units->"Significance Probability"]
]],
If[Length[OptionValue[\[ScriptCapitalH]2]]==0,"",
tableout[Quiet[Map[{#["TestDataTable",All],#["ShortTestConclusion"]}&,Optio=
nValue[\[ScriptCapitalH]2]]],
OptionValue[desc]<>": Test \[ScriptCapitalH]2 that the variance of the=
original & recreated
variables are the same (ie that " <> ToString[k] <>" factors is
sufficient)",
varnames,{"Statistic","Conclusion"},units->"Variables"]
],
tableout[Correlation[\[CapitalLambda]],OptionValue[desc]<>": Correlation=
of Factors",facnames,facnames,units->"R"],
tableout[Covariance[\[CapitalLambda]],OptionValue[desc]<>": Covariance of=
Factors",facnames,facnames,units->"Cov"],
tableout[\[CapitalLambda]5,OptionValue[desc]<>": Factor Loadings Squares:=
Variance in each variable explained
by factor & Relative
Importance",Flatten[{varnames,{"Total","Ratio"}}],Flatten[{facnames,{"Commu=
nality"},rinames}]],
tableout[Transpose[{\[CapitalLambda]2,Rest[FoldList[Plus,0,\[CapitalLambda]=
2]]}],OptionValue[desc]<>":Factor
Sum of the Squares",facnames,{"Variance
Proportion","Cumulative
Variance"}],
If[Length[OptionValue[fs]]==0,"",
tableout[OptionValue[fs],OptionValue[desc]<>": Factor
Scores",obsnames,facnames]
],
"",
If[k>=3,ListPlot3D[Take[\[CapitalLambda],All,3],ImageSize->Medium,PlotLab=
el->OptionValue[desc]<>":
Plot of First 3 Factors for all Variables"],""],
"",
If[k>=2,ListPlot[Take[\[CapitalLambda],All,2],ImageSize->Medium,PlotLabel=
->OptionValue[desc]<>":
Plot of First 2 Factors for all Variables"],""]
}]]
(*Ratio of factor eigenvalues - explanatory importance of the factors with=
respect to the total sample*)
(*{\[CapitalLambda]u,\[CapitalLambda]d,\[CapitalLambda]v}=SingularValueDe=
composition[\[CapitalLambda],k];
tableout[Sqrt[Diagonal[\[CapitalLambda]d]]/Total[\[CapitalLambda]d,2],0.001=
,"Ratio
of factor eigenvalues",facnames,"Proportion",""]*)
];
varimaxTP[\[CapitalLambda]_,OptionsPattern[{NormalizeVars->False}]]:=Modu=
le[{iter,itermax,\[Epsilon],amat,amatd,p,k,targetbasis,rotm,varnow,notconve=
rged,phimax,subrot,varold,optamat,scalefactor},
(*varimaxTP(component loadings matrix) gives an unnormalized Varimax
optimal rotation of matrix AMAT using the standard algorithm
Kaiser(Psychometrika 1958) modified by Trevor Park in April 2002 from an=
original file by J.O.Ramsay. See
http://www.stat.ufl.edu/~tpark/Research/varimaxTP.m*)
itermax=100;(*Maximum number of complete passes through column pairs*)
\[Epsilon]=10^-7;(*Stopping tolerance based on relative reduction in=
varimax*)
If[Length[Dimensions[\[CapitalLambda]]]!=2, Print["Error: AMAT must be=
two-dimensional"];Return[\[CapitalLambda]]];
{p,k}=Dimensions[\[CapitalLambda]];
If[OptionValue[NormalizeVars],
amat=Map[Normalize,\[CapitalLambda]],(*Normalize the variables (rows) of=
\[CapitalLambda]*)
amat=\[CapitalLambda]];
rotm=IdentityMatrix[k];
targetbasis=IdentityMatrix[p];
varnow=Total[Variance[amat^2]];
notconverged=True;
iter=0;
While[notconverged && iter<=itermax,
For[j=1,j<=(k-1),j++,
For[l=(j+1),l<=k,l++,
(*Calculate optimal 2-D planar rotation angle for columns j,l*)
phimax=Arg[p
Total[Thread[Complex[amat[[All,j]],amat[[All,l]]]]^4]-Total[Thread[Complex[=
amat[[All,j]],amat[[All,l]]]]^2]^2]/4;
subrot={{Cos[phimax],-Sin[phimax]},{Sin[phimax],Cos[phimax]}};(*orthonorm=
al
matrix*)
amat[[All,{j,l}]]=amat[[All,{j,l}]].subrot;
rotm[[All,{j,l}]]=rotm[[All,{j,l}]].subrot;
];
varold=varnow;
varnow=Total[Variance[amat^2]];
If[varnow==0,Break[]];
notconverged=(varnow-varold)/varnow>\[Epsilon];
iter++
];
If[OptionValue[NormalizeVars],(*Rescale normalize the variables (rows) of=
\[CapitalLambda]*)
scalefactor=Sqrt[Total[\[CapitalLambda]^2,{2}]]/Sqrt[Total[amat^2,{2}]];
Print[Dimensions[scalefactor]];
amat=Map[amat[[#]]*scalefactor[[#]]&,Range[p]]];
If[iter>=itermax,Print["Maximum number of iterations reached in
function!"]];
optamat=targetbasis.amat;
Return[optamat]
];
];
End[]
EndPackage[]
UTS CRICOS Provider Code: 00099F
DISCLAIMER: This email message and any accompanying attachments may contain=
confidential information.
If you are not the intended recipient, do not read, use, disseminate, distr=
ibute or copy this message or
attachments. If you have received this message in error, please notify the =
sender immediately and delete
this message. Any views expressed in this message are those of the individu=
al sender, except where the
sender expressly, and with authority, states them to be the views of the Un=
iversity of Technology Sydney.
Before opening any attachments, please check them for viruses and defects.
Think. Green. Do.
Please consider the environment before printing this email.
- References:
- Compare matrices using Chi square for Exploratory Factor Analysis
- From: Stuart Nettleton <Stuart.Nettleton@uts.edu.au>
- Compare matrices using Chi square for Exploratory Factor Analysis