*! version 1.3.0 Philippe VAN KERM, February 1998
*! Multiple Correspondence Analysis
*! Syntax: mca <varlist>  [if <exp>] [in <range>] [weight] [,d(#) q(#) notrans]


capture program drop mca
program define mca


version 5.0
set more on

loc varlist "req ex min(2)"
loc in "opt"
loc if "opt"
loc weight "aweight fweight"
loc options "d(integer 0) q(real 0) Notrans"
parse "`*'"
if `d'<0 | `d'==1 {
	di in red "invalid number of dimensions"
	exit 198
	}
if `q'<0 | `q'>1 {
	di in gr "Quality threshold mis-specified : option q(`q') rejected"
	loc q = 0
	}
loc weight "[`weight'`exp']"
parse "`varlist'", parse (" ")
loc n1 1
while "`1'" ~= "" {
	loc var`n1' `1'
	capture confirm string variable `var`n1''
	if _rc==0 {
		di in red "variables in varlist must be numeric"
	exit 109
	}
	mac shift
	loc n1 = `n1' + 1
	}
loc nbv = `n1' - 1

tempname  B N IN JN tmp1 tmp2 RS CS RM CM P DR DC DR12 DC12 tmp3 A PrinIn /*
       */ TotIn U T V OT UU VV DL1 DL2 DL1C DL2C RI CI SQC1 SQC2 F FP GP n /*
       */ PrinInp tmp

/*** Burt matrix ***/

preserve

loc n1 1
while `n1'<=`nbv' {
	qui drop if `var`n1''==.
	loc n1 = `n1' + 1
	}

loc n1 1
while `n1' <= `nbv' {
	loc n2 = 1
	while `n2' <= `nbv' {
		qui tab `var`n1'' `var`n2'' `weight' `if' `in' , /*
		    */ matcell(N`n1'to`n2') matrow(IN`n1'to`n2') /*
		    */ matcol(JN`n1'to`n2')
		if _result(1)==0 {
			di in blue "no observations"
			exit
			}
		loc n2 = `n2' + 1
		}
	loc n1 = `n1' + 1
	}
loc n1 1
while `n1' <= `nbv' {
	loc n2 2
 	mat B`n1' = N`n1'to1
	mat drop N`n1'to1
	while `n2' <= `nbv' {
		mat B`n1' =  B`n1',N`n1'to`n2'
		mat drop N`n1'to`n2'
		loc n2 = `n2' + 1
		}
	loc n1 = `n1' + 1
	}
loc n1 2
mat `B' = B1
mat drop B1
while `n1' <= `nbv' {
	mat `B' =  `B'\B`n1'
	mat drop B`n1'
	loc n1 = `n1' + 1
	}

/*** Correspondence analysis on the Burt matrix ***/
mat `N' = `B'
loc i = rowsof(`N')
loc j = colsof(`N')
mat `RS' = J(`i',1,0)
mat `CS' = J(1,`j',0)
loc n1 1
loc nobs 0
while `n1' <= `i' {
	loc n2 1
	loc sum 0
	while `n2' <= `j' {
		loc sum = `sum' + `N'[`n1',`n2']
		loc n2 = `n2' + 1
		}
	mat `RS'[`n1',1] = `sum'
	loc nobs = `nobs' + `sum'
	loc n1 = `n1' + 1
	}
scalar n = `nobs'
loc n1 1
while `n1' <= `j' {
	loc n2 1
	loc sum 0
	while `n2' <= `i' {
		loc sum = `sum' + `N'[`n2',`n1']
		loc n2 = `n2' + 1
		}
	mat `CS'[1,`n1'] = `sum'
	loc n1 = `n1' + 1
	}
scalar overn = 1/scalar(n)
mat `RM' = overn*`RS'
mat `CM' = overn*`CS'
mat `P' = overn*`N'
mat `DR' = diag(`RM')
mat `DC' = diag(`CM')

/***     Computation of the matrix of standardized residuals     ***/

mat `DR12'=J(`i',`i',0)
loc n1 1
while `n1'<=`i' {
	mat `DR12'[`n1',`n1'] = (1/(sqrt(`DR'[`n1',`n1'])))
	loc n1 = `n1'+1
	}
mat `DC12'=J(`j',`j',0)
loc n1 1
while `n1'<=`j' {
	mat `DC12'[`n1',`n1'] = (1/(sqrt(`DC'[`n1',`n1'])))
	loc n1 = `n1'+1
	}
mat `tmp1' = `RM'*`CM'
mat `tmp1' = `P'-`tmp1'
mat `tmp1' = `DR12'*`tmp1'
mat `A' = `tmp1'*`DC12'

/***     Singular Value Decomposition and Computations of Inertia   ***/

/* Singular value decomposition */

mat svd `U' `T' `V' = `A'

/* Ordering of U T and V */

mat `OT' = `T'
mat `UU' = `U'
mat `VV' = `V'
loc n1 1
while `n1' <= `j' {
	loc n2 1
	loc max =  -1
	loc posmax 0
	while `n2' <= `j' {
		if `T'[1,`n2'] > `max' {
			loc max = `T'[1,`n2']
			loc posmax = `n2'
			}
		loc n2 = `n2' + 1
		}
	mat `T'[1,`posmax'] = 0 - 1
	mat `OT'[1,`n1'] = `max'
	loc n3 1
	while `n3'<=`i' {
		mat `UU'[`n3',`n1'] = `U'[`n3',`posmax']
		loc n3 = `n3' + 1
		}
	loc n3 1
	while `n3'<=`j' {
		mat `VV'[`n3',`n1'] = `V'[`n3',`posmax']
		loc n3 = `n3' + 1
		}
	loc n1 = `n1' + 1
	}
mat `U' = `UU'
mat `V' = `VV'
mat `T' = diag(`OT')

loc K = `j' - `nbv'
mat `PrinIn' = J(3,`K',0)
loc n1 1
while `n1' <= `K' {
	mat `PrinIn'[1,`n1'] = (`T'[`n1',`n1'])^2
	loc n1 = `n1' + 1
	}

if "`notrans'" == "" {
	loc n1 1
	while `n1' <= `K' & sqrt(`PrinIn'[1,`n1'])>1/`nbv' {
		mat `PrinIn'[1,`n1'] = /*
                 */ ((`nbv'/(`nbv'-1))^2) * (((`T'[`n1',`n1'])-(1/`nbv'))^2)
		loc n1 = `n1' + 1
		}
	loc newK = `n1'-1
	while `n1' <= `K' {
		mat `PrinIn'[1,`n1'] = 0
		loc n1 = `n1' + 1
		}
	loc K = `newK'
	}
scalar TotIn = 0
loc n1 1
while `n1' <= `K' {
	scalar TotIn = scalar(TotIn) + `PrinIn'[1,`n1']
	loc n1 = `n1' + 1
	}
loc n1 1
loc sum 0
while `n1' <= `K' {
	mat `PrinIn'[2,`n1'] = `PrinIn'[1,`n1']/scalar(TotIn)
	loc sum  = `sum' + `PrinIn'[2,`n1']
	mat `PrinIn'[3,`n1'] = `sum'
	loc n1 = `n1' + 1
	}

loc n1 1
while `n1' <= `K' {
	mat `T'[`n1',`n1'] = sqrt(`PrinIn'[1,`n1'])
	loc n1 = `n1' + 1
	}

/* Row principal coordinates in F and decompositions of inertia by rows  */

mat `tmp1' = `DR12'*`U'
mat `tmp2' = `tmp1'*`T'
mat `F' = J(`i',`K',0)
mat `DL1' = J(`i',`K',0)
mat `DL1C' = J(`i',`K',0)
loc n1 1
while `n1' <= `i' {
	loc n2 1
	while `n2' <= `K' {
		mat `F'[`n1',`n2'] = `tmp2'[`n1',`n2']
		mat `DL1'[`n1',`n2'] = ((`F'[`n1',`n2'])^2)*(`RM'[`n1',1])
		mat `DL1C'[`n1',`n2'] = `DL1'[`n1',`n2'] / `PrinIn'[1,`n2']
	 	loc n2 = `n2' + 1
		}
	loc n1 = `n1' + 1
	}

mat `RI' = J(`i',1,0)
loc n1 1
while `n1' <= `i' {
	loc n2 1
 	loc sum 0
	while `n2'<=`K' {
		loc sum = `sum' + `DL1'[`n1',`n2']
		loc n2 = `n2' + 1
		}
	mat `RI'[`n1',1] = `sum'
	loc n1 = `n1' + 1
	}
mat `SQC1'=J(`i',`K',0)
loc n1 1
while `n1' <=`i' {
	loc n2 1
	while `n2' <= `K' {
		mat `SQC1'[`n1',`n2'] = `DL1'[`n1',`n2'] / `RI'[`n1',1]
		loc n2 = `n2' + 1
		}
	loc n1 = `n1' + 1
	}

/**************     Output Display     **************/

title

/*** Table display ***/

loc n1 1
loc n2 1
loc rname1 " "
while `n1' <= `nbv' {
	loc le = rowsof(IN`n1'to`n1')
	loc n3 1
	while `n3' <= `le' {
		loc nam = IN`n1'to`n1'[`n3',1]
		loc ncar = length("`var`n1''")
		if `ncar'>5 {loc ncar 5}
		loc name = substr("`var`n1''",1,`ncar') + "_`nam'"
		if  length("`name'")>8 {loc name = substr("`name'",1,8)}
                loc rname`n2' = "`rname`n2''" +  "`name'" + " "
		if length("`rname`n2''")>=50 {loc n2 = `n2' + 1}
		loc n3 = `n3' + 1
		}
	loc n1 =`n1' + 1
	}

if (`d'==0) | (`d'>`K') {loc dg = `K'}
else {loc dg = `d'}

loc n1 1
loc n2 1
loc xname1 " "
while `n1' <=`dg' {
	loc nam = "Dim`n1'"
	loc xname`n2' = "`xname`n2''" + "`nam'" + " "
	if length("`xname`n2''")>=50 {loc n2 = `n2' + 1}
	loc n1 =`n1' + 1
	}

di
di in green " Total Inertia : " in yellow %10.3f scalar(TotIn)
di

mat `PrinInp' = `PrinIn'[.,1..`dg']
mat `PrinInp' = `PrinInp''
mat colnames `PrinInp' = Inertia Share Cumul
mat rownames `PrinInp' = `xname1'`xname2'`xname3'`xname4'`xname5' /*
                    */ `xname6'`xname7'`xname8'`xname9'`xname10'
di in green " Principal Inertia Components :"
mat list `PrinInp', noheader f(%10.3f)


mat `tmp' = `F'[.,1..`dg']
mat `tmp' = `RI',`tmp'
mat `FP' = `RM',`tmp'
mat rownames `FP' =`rname1'`rname2'`rname3'`rname4'`rname5' /*
              */ `rname6'`rname7'`rname8'`rname9'`rname10'
mat colnames `FP' = Mass  Inertia `xname1'`xname2'`xname3'`xname4'`xname5' /*
               */ `xname6'`xname7'`xname8'`xname9'`xname10'
di
di in green " Coordinates : "
mat list `FP' , noheader f(%10.3f)

mat `tmp' = `DL1C'[.,1..`dg']
mat rownames `tmp' = `rname1'`rname2'`rname3'`rname4'`rname5' /*
                */ `rname6'`rname7'`rname8'`rname9'`rname10'
mat colnames `tmp' = `xname1'`xname2'`xname3'`xname4'`xname5' /*
                */ `xname6'`xname7'`xname8'`xname9'`xname10'
di
di in green " Explained inertia of axes  :"
mat list `tmp' , noheader f(%10.4f)

mat `tmp' = `SQC1'[.,1..`dg']
mat rownames `tmp' = `rname1'`rname2'`rname3'`rname4'`rname5' /*
                */ `rname6'`rname7'`rname8'`rname9'`rname10'
mat colnames `tmp' = `xname1'`xname2'`xname3'`xname4'`xname5' /*
                */ `xname6'`xname7'`xname8'`xname9'`xname10'
di
di in green " Contributions of principal axes :"
mat list `tmp' , noheader f(%10.4f)

/*** Graphical display ***/

if `d' ~= 0 {
	if `K' < `d' {
di in red "Not enough principal dimensions for the representations requested"
        }
	else {
		mat colnames `F' = `xname1'`xname2'`xname3'`xname4'`xname5' /*
                              */ `xname6'`xname7'`xname8'`xname9'`xname10'
		mat F = `F'[.,1..`dg']
		svmat F ,names(matcol)
		qui gen str20 Flab=" "
		loc n1 1
		loc n3 1
		while `n1' <= `nbv' {
			loc le = rowsof(IN`n1'to`n1')
			loc n2 1
			while `n2' <= `le' {
				loc nam = IN`n1'to`n1'[`n2',1]
				loc ncar = length("`var`n1''")
				if `ncar'>5 {loc ncar 5}
qui replace Flab=substr(substr("`var`n1''",1,`ncar') + "_`nam'", 1,8) in `n3'
				loc n2 = `n2' + 1
				loc n3 = `n3' + 1
				}
			loc n1 =`n1' + 1
			}
		if `q'~=0 {
			loc n1 1
			while `n1'<=`i' {
				loc n2 1
				loc qual 0
				while `n2' <=`dg' {
					loc qual = `qual' + `SQC1'[`n1',`n2']
					loc n2 = `n2' + 1
					}
                             if `qual' < `q' {qui replace Flab="." in `n1'}
				loc n1 =`n1' + 1
				}
			}
		loc n1 = 1
		while `n1' <= `d'-1 {
			loc n2 = `n1' + 1
			while `n2'<=`d' {
                               graph FDim`n2' FDim`n1', s([Flab]) /*
                                  */ b1("Multiple Correspondence Analysis") /*
                                  */ b2("Axis `n1'") l1(" ") /*
                                  */ l2("Axis `n2'") xla yla xline(0) /*
                                  */ yline(0)  border
				loc n2 = `n2' + 1
				more
				}
			loc n1 = `n1' + 1
			}

		}
	}

restore

end

capture program drop title
program define title
	di in green _d(78) "-"
	di
	di in green _col(20) "MULTIPLE CORRESPONDENCE ANALYSIS"
	di
	di in green _d(78) "-"
	di
end