* promax.sps------------------------------------------------------------.
* 因子分析（主因子解→プロマックス回転）マクロ for spss. by k.hori(hori@ec.kagawa-u.ac.)
*   03/ 3/26 Wed.
* var=変数リスト(必須)/nfact=因子数（既定値３）/k=プロマックス累乗係数(既定値３)
* /iterate=最大反復数(既定値30)/tol=精度(既定値0.0001).
*-------------------------------------------------------------.
define hfactor (var=!charend('/') 
/nfact=!default(3)!charend('/')
/k=!default(3)!charend('/')
/iterate=!default(100)!charend('/')
/tol=!default(0.0001)!charend('/')
).
preserve.
set length=none printback = off   mxloops = 90000. 
matrix.

*因子数.
compute nfactors=!nfact.
*promax k.
compute ppower=!k.
*反復数.
compute itercfa=!iterate.
*精度.
compute tolprn=!tol.

get data / file=* /missing=omit/ name=varname1 /variables = !var.

compute cases = nrow(data). 
compute nvars = ncol(data).
compute varname1={varname1,'2乗和'}.
print {cases,nfactors,ppower,itercfa}/title = '階層因子分析(主因子解→promax 回転)'
  /clabels='N','因子数','k','最大反復数'/format=f5.0.
print tolprn/title="許容限度(ε)"/format=f8.7.
compute varnames=varname1.

*主因子解.
compute extract=1.
*プロマックス回転.
compute rotate=1.
compute core=nvars.

compute num2={' 1',' 2',' 3',' 4',' 5',' 6',' 7',' 8',' 9',
'10','11','12','13','14','15','16','17','18','19',
'20','21','22','23','24','25','26','27','28','29',
'30','31','32','33','34','35','36','37','38','39',
'40','41','42','43','44','45','46','47','48','49',
'50','51','52','53','54','55','56','57','58','59',
'60','61','62','63','64','65','66','67','68','69',
'70','71','72','73','74','75','76','77','78','79',
'80','81','82','83','84','85','86','87','88','89',
'90','91','92','93','94','95','96','97','98','99','100'
}.
compute fs={"F 1","F 2","F 3","F 4","F 5","F 6","F 7","F 8","F 9","F10","F11","F12","F13","F14","F15"}.
compute hos={"HO 1","HO 2","HO 3","HO 4","HO 5","HO 6","HO 7","HO 8"}.

* 相関行列.
* correlation matrix.
compute zscores = make(nrow(data), ncol(data), -999).
compute n   = nrow(data).
compute nm1 = n - 1.
loop #a = 1 to ncol(data).
  compute mean = csum(data(:,#a)) / n.
  compute sd = sqrt ( cssq((data(:,#a)- mean)) / nm1  ).
  compute zscores(:,#a) = (data(:,#a) - mean) / sd.
end loop.
compute rdata = (1 / (nrow(data)-1) ) * ( t(zscores) * zscores ) .

compute rcore =  rdata(1:core,1:core) .
compute r = rcore.

*階層因子分析のためのループ.
loop #ll= 1 to 5.
do if (#ll>1).
   print (#ll-1)/title='高次因子'.
end if.
compute evals1 = eval(r).
print t(evals1)/title='対角１の相関行列の固有値'/format=f8.4.
* Factor Extraction    >>>>>>.
* 主因子解   extract=1; 最大反復数 itercfa; 相関行列 r; =>負荷量 lding.
* CFA / PAF  -- from Bernstein p 189   --  smc = from Bernstein p 104.
do if ( extract = 1).
compute rcfa = r.
compute smc = 1 - (1 &/ diag(inv(rcfa)) ).
loop #a = 1 to itercfa.
call setdiag(rcfa,smc).
call eigen(rcfa,eigvect,eigval).
compute eigval = mdiag(eigval).
compute lding=eigvect(:,1:nfactors)*sqrt(eigval(1:nfactors,1:nfactors)).
compute communal = rssq(lding).
do if ( mmax(abs(communal-smc))< tolprn ).
break.
else.
compute smc=communal.
end if.
end loop.
end if.
do if #a< itercfa.
   print #a/title='反復数'.
else.
   print #a/title='収束しませんでした：反復数'.
end if.  

** 固有ベクトルの正負の入れ替え.
compute fugo=csum(lding).
compute fugo=(fugo>=0)+((fugo<0)*(-1)).
loop #b = 1 to nfactors.
compute lding(:,#b)=lding(:,#b)*fugo(#b).
end loop.

compute names=num2.
compute names((nfactors+1))='共通性'.
print {lding,smc}/rnames=varnames/title= '主因子法　因子行列'/format= f8.3/cnames=names.


* 主成分分析終了.
* バリマックス回転→プロマックス回転  rotate=1;.
* Promax rotation -- Marcus, 1993.
do if ( rotate = 1 and nfactors>1).
* varimax rotation -- Marcus, 1993, in Reyment & Joreskog, 
  Applied Factor Analysis in the Natural Sciences, CUP.
  compute b=lding.
  compute n = nrow(lding).
  compute nf = ncol(lding).
  compute hjsq=diag(lding*t(lding))   .
  compute hj=sqrt(hjsq).
  compute bh=lding &/ (hj*(make(1,nf,1)))  .
  compute Vtemp=n*rsum(csum(bh&**4))-rsum(csum(bh&**2)&**2)  .
  compute V0=Vtemp.
  loop #it=1 to 20 . 
    loop #i =1 to nf-1 .
      compute jl=#i+1 .
      loop #j=jl to nf.
        compute xj=lding(:,#i) &/ hj.   
        compute yj=lding(:,#j) &/ hj.   
        compute uj=xj &* xj-yj &* yj.
        compute vj=2*xj &* yj.
        compute A=csum(uj).
        compute bigB=csum(vj).
        compute C=t(uj) * uj - t(vj) * vj.
        compute D=2 * t(uj) * vj.
        compute num=D-2*A*bigB/n.
        compute den=C-(A**2-bigB**2)/n.
        compute tan4p=num/den.
        compute tan4p = num / den .
        compute artanp = ( artan(abs(tan4p)) ) / 4.
        do if (num>0 and den>0).
           compute phi = artanp.
        else if (num<0 and den<0).
             compute phi =  -1 * (.7854 - artanp)  .
           else if (num<0 and den>0).
               compute phi =  -1 * artanp.
             else if (num>0 and den<0).
                compute phi =  .7854 - artanp .
       end if.
       compute angle = phi * 180 / 3.1415926  .
 *    do if (  abs(phi)>.00001  ).
      do if (  abs(phi)>tolprn  ).
         compute bigXj=cos(phi)*xj+sin(phi)*yj.
         compute bigYj=-sin(phi)*xj+cos(phi)*yj.
         compute bj1=bigXj &* hj.
         compute bj2=bigYj &* hj.
         compute b(:,#i)=bj1.
         compute b(:,#j)=bj2.
         compute lding(:,#i)=b(:,#i).
         compute lding(:,#j)=b(:,#j).
       end if.
      end loop.
    end loop.
    compute lding=b.
    compute bh=lding &/ (hj*(make(1,nf,1)))  .
    compute Vtemp=n*rsum(csum(bh&**4))-rssq(cssq(bh))  .
    compute V=Vtemp.
    do if ( abs(V-V0)< tolprn ).
      break.
    else.
      compute V0=V.
    end if.
  end loop.
  print {lding;cssq(lding)}/rnames=varnames/title= 'バリマックス回転　因子負荷量'/format= f8.3/cnames=num2.
* プロマックス回転 rotate=1; ppower 累乗.
  compute B = lding.
*目標行列がおかしいようなので修正する、.
  compute bstar=make(n,nf,0).
*共通性.
  compute comm=rssq(B).
  compute comm=sqrt(comm).
  loop #i=1 to n.
    compute bstar(#i,:)=B(#i,:)&/comm(#i).
  end loop. 
  compute w=cmax(abs(bstar)).
  loop #i =1 to n.
    compute ww=Bstar(#i,:)&/w.
    compute sgn=((ww>0)-(ww<0)).
    compute bstar(#i,:)=(abs(ww)&**ppower)&*sgn.
  end loop.
 * print Bstar/rnames=varnames/title= '目標行列'/format= f8.3/cnames=num2.
  compute Tr=inv(t(B) *B)*t(B)*Bstar.
 * print tr/title="プロクラステス変換行列"/format=f8.3.
  compute Tr=Tr*sqrt(inv(mdiag(diag(t(Tr)*Tr)))).
  compute Tpp=inv(Tr).
  compute Tpp=inv(sqrt(mdiag(diag(Tpp*t(Tpp)))))*Tpp.
*参考構造.
  compute Sr=B*Tr.
*print Sr/rnames=varnames/cnames=num2/title='参考構造'/format= f8.3.
*因子間相関.
  compute Phip=Tpp*t(Tpp).
*因子構造行列.
  compute Sp=B*t(Tpp).
*因子パタン行列.
  compute Pp=B*inv(Tpp).
  print {Pp;cssq(Pp)}/rnames=varnames/cnames=num2/title='プロマックス法　因子パタン行列'/format= f8.3.
  print Phip/rnames=num2/cnames=num2/title='因子間相関'/format= f8.3.
  print Sp/rnames=varnames/cnames=num2/title='因子構造行列'/format= f8.3.

*参考パタン行列.
  compute Pr=B*inv(t(Tr)).
 * print Pr/rnames=varnames/cnames=num2/title='参考パタン行列'/format= f8.3.
*参考軸間相関.
  compute Phir= t(Tr)*Tr.
 * print Phir/rnames=num2/cnames=num2/title='参考軸間相関'/format= f8.3.
*高次因子分析のための相関行列.
* プロマックス回転終了.
end if.
do if (#ll=1).
   do if (nfactors<2).
       break.
    end if.
 *  print /title='高次因子'.
   compute lding1=lding.
   compute varnames=fs.
   compute varnames((nfactors+1))='2乗和'.
   compute nfactors=1.
   compute r=Phip.  
else.
   compute rawnames={hos(1:nfactors),fs}.
   compute d=mdiag(sqrt(1-smc)).
   compute dd={lding,d}.
   compute ldhigh=lding1*(dd).
   compute com=rssq(ldhigh).
   print {ldhigh;cssq(ldhigh)}/format=f8.3/rnames=varname1/cnames=rawnames/title='階層因子分析 因子パタン'.
   compute nfactors=nfactors+1.
   do if (evals1(#ll)<=1).
       break.
    end if.

end if.

end loop. 
/* #ll 階層因子分析のためのループ終了 */.
end matrix.
restore.
!enddefine.
hfactor var=v1 to v24/nfact=4.