/*
/ Program   : Exploratory factor analysis
/ Version   : 1.0
/ Editor(s) : Nicole M Lindner
/ Date      : May 6, 2009
/ Contact   : nml5d@email.virginia.edu
/===========================================================================
/ Purpose   : Factor Analysis
/	*	Output PROC CORR values to new SAS dataset
/	* 	Conduct an exploratory factor analysis, per Russell (2002)'s 
/			recommendations
/	*	As such, calculate the predicted eigen values with bootstrap (?)
/			sampling of randomly-distributed data of the same dimensions
/			as my sample
/	THIS IS USELESS WITHOUT THE FOLLOWING REFERENCE:
/	Russell, D. W. (2002).  In search of underlying dimensions: The use 
/	  (and abuse) of factor analysis in PSPB. Personality and Social 
/	  Psychology Bulletin, 28, 1629-1646.
/	Available at: http://psp.sagepub.com/cgi/reprint/28/12/1629
***************************************************************************
/ AMENDMENT HISTORY:
/ init --date--- mod-id ----------------------description------------------
/ 1.0 May5,2009 Documentation:
/					Adds documentation before adding script to Dataverse
/
/==========================================================================
/ This is public domain software. No guarantee as to suitability or 
/ accuracy is given or implied. 
/ User uses this code entirely at their own risk.
/=========================================================================*/


/* ----------- Outputting PROC CORR values to new SAS dataset ----------- */
	* Outputs the Pearson correlations between my variables of interest in 
	   a dataset that's structured for use with PROC FACTOR.;
	* Note: The NOMISS means that only observations with complete data 
	   will be included;

DATA CRS1;SET CR.FinalCR;
LABEL q01=	"1  Should allow Actions";
LABEL q02=	"2R Should prohibit Actions";
LABEL q03=	"3  Const protects actions";
LABEL q04=	"4R Const prohibits actions";
LABEL q05=	"5R Bad idea to express this way";
LABEL q06=	"6R Shouldnt have expressed opinion at all";
LABEL q07=	"7  Has the right to express his opinion";
LABEL q08=	"8R Pass law to prevent speech like this";
LABEL q09=	"9R Const should restrict";
LABEL q10=	"10 Const should protect";
LABEL q11=	"11R People may be harmed";
LABEL q12=	"12R oppose his opinion";
LABEL q13=	"13 agree w/ his opinion";
run;

	* The 11 items assessing willingness to protect the speech;
PROC CORR DATA=CrS1 OUTP=IniFactor NOMISS;	VAR Q01-Q11;run; 

	* Initial checking to make sure output worked;
PROC PRINT DATA=IniFactor;RUN; ; 

/* --------- Exploratory factor analysis, per Russell (2002)    --------- */
PROC FACTOR DATA=IniFactor (TYPE=CORR)
                METHOD=PRINCIPAL	PRIORS=SMC	MINEIGEN=1
                CONV= 0.001	MAXITER=200
                MSA	SIMPLE	SCREE	REORDER	HEYWOOD;
RUN;
/* --------- Generating Average eigen values for the sample size--------- */
/*	This is adapted from
	the script in Russell (2002), Footnote 6. As an R novice, the script
	in Russell was difficult to follow. Russell recommended that a full 
	script be downloaded from the R-Project's source files, but I still 
	can't locate it. Per Russell's citation, I consulted Reise, Waller,
	& Comfrey (2000), but had trouble with their R code. Thus, I've 
	simplified its use for SAS users. 
	Copy and paste the codes after the $x into R. 
	Replace all references to:
	 * 11 with your number of variables, 
	 * 1706 with your N, and 
	 * the real.eig vector with your eigenvalues (from SAS or another 
		statistical program).
	This then plots your data's eigen values against the  bootstrapped
	average null eigen values for a datset of your dimensions.
$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
###	This will also plot your eigen values as a thick blue line & solid dots,
#	against the empty red dots of the average null eigen values for a
#	sample of your size.
random.eig <- matrix(0,nrow =100, ncol = 11) 
	for(i in 1:100) 	{
		random.data <- matrix(rnorm(1706 * 11), nrow = 1706, ncol = 11) 
		random.eig[i, ] <- eigen(cor(random.data)) $values
		}

Nfactors<-c(1:11)
Nfactors
average.eig <- apply(random.eig,2,mean)
average.eig
random.eig

###	Eigen-values from PROC FACTOR in SAS, using communalities as the 
	initial estimates. I'm simply reading in my data's eigen values
	that SAS reported from PROC FACTOR. Change these to match yours.
real.eig <-c(4.78653937,0.63354917,0.25127438,0.18678884,0.08918386,-.00325701,-.09540982,-.12071150,-.13477573,-.18154008,-.19962807)
#	same here, change 1:11 in the plot and lines statements with your 
#	number of factors, change the limits of the y axis (-2 to 5) to fit 
#	your ranges

###	Plots them against each other. Change the c(1:X) in the plot and 
#	lines statements to fit the number of factors
plot(1:11, real.eig, type="o", col="blue",cex.lab= 1.3,
	xlab="Number of Factors",
	ylab="Eigen Values",
	cex=1.25, pch=19,lwd=3,lty=1)
axis(1, at=1:11,labels=FALSE,tick=TRUE)
lines(1:11,average.eig,type="b",pch=21,col="red",lwd=.5,lty=2)
abline(h=1,untf=FALSE)

legend("topright", c("'Real' Eigen Values from Data","Random Eigens from Parallel Analysis","1=explains more variance (like F >= 1)"), 
	col = c("blue","red","black"),cex=1.3,
       lty = c(1,2), pch = c(19,21,NA_integer_),
       merge = TRUE)
#	lty = line type and width for lines in legend
#	pch = the plotting symbols appearing in the legend 
*************/