*! 1.0.1 Stephen Soldz 11/12/2002 program define minap, rclass version 7 ************************************************************************** * Calulates Velicer(1976) minimum average partial correlation ** * estimate of number of principal components to extract. ** * Author: Stephen Soldz, Boston Graduate School of Psychoanalysis, ** * 11/8/2002 ** * Reference: Velicer, W.F. (1976). Determining the number of components ** * from the matrix of partial correlations. Psychometrka, 41, 321-327. ** ************************************************************************** capture syntax varlist(min=3 numeric) [if] [in] if _rc { capture syntax , corr(str) } if _rc { display as err "specify either varlist of at least 3 variables or corr()" exit 198 } tempname Cov N R Rstar X v sqrv A p pminus1 Cstar i j k fm f0 fmvec tempname numeigen index prior if "`varlist'" ~= "" { qui mat accum `Cov' = `varlist' `if' `in', dev noconst local N = r(N) mat `Cov' = `Cov' / (`N' - 1) mat `R' = corr(`Cov') } else if "`corr'" ~= "" { mat `R' = `corr' } local p = colsof(`R') scalar `f0' = 0 forval i = 2 / `p' { local j = `i' - 1 forval k = 1 / `j' { scalar `f0' = `f0' + `R'[`i', `k']^2 } } scalar `f0' = 2 * `f0' / (`p' * (`p' - 1)) dis as text _ne "{title:Minimum Average Partial Correlation for Number of Principal Components}" dis as text _ne "NOTE: Pick number of components (m) at which fm is minimum." dis as text " If f1 > f0 (average intervariable correlation) " dis as text " then no components should be extracted." dis as text _ne _col(5) "m = 0" _col(15) "f0 = " _col(25) as result `f0' _ne * Initialize holder of fm values mat `fmvec' = J(1, `p' - 1, 0) mat symeigen `X' `v' = `R' * get number of eigenvalues > 1 local numeigen = 0 forval i = 1 / `p' { local numeigen = `numeigen' + (`v'[1, `i'] > 1.0) } * Rescale `v' `X' so square of columns sum to eigenvalues mat `sqrv' = J(`p', `p', 0) forval i = 1 / `p' { mat `sqrv'[`i', `i'] = sqrt(`v'[1, `i']) } mat `X' = `X' * `sqrv' local pminus1 = `p' - 1 forval m = 1 / `pminus1' { ** Main loop for calculating mat `A' = `X'[1..`p', 1..`m'] mat `Cstar' = `R' - `A' * `A'' mat `Rstar' = corr(`Cstar') scalar `fm' = 0 forval i = 2 / `p' { local j = `i' - 1 forval k = 1 / `j' { scalar `fm' = `fm' + (`Rstar'[`i', `k']^2) } } scalar `fm' = 2 * `fm' / (`p' * (`p' - 1)) mat `fmvec'[1, `m'] = `fm' dis _col(5) as text "m = " `m' _col(15) "f`m' = " /* */ _col(25) as result `fm' } local index = 0 scalar `prior' = `f0' forval i= 1 / `pminus1' { if `fmvec'[1, `i'] < `prior' { scalar `prior' = `fmvec'[1, `i'] local index = `i' } } local s = cond(`index' != 1, "s", "") dis as text _ne "{p}{cmd:minap} procedure suggests that {res:`index'} " dis as text "principal component`s' should be extracted" local s = cond(`numeigen' != 1, "s", "") dis as text _ne "{p}For comparison, the Kaiser eigenvalue > 1 rule suggests extracting" dis as text "{res:`numeigen'} principal component`s'" return scalar minap =`index' return scalar Kaiser = `numeigen' return matrix Cormat `R' return matrix EigenVal `v' return matrix EigenVec `X' return matrix FmVec `fmvec' end