topic Re: Factor analysis using maximum likelihood estimation in SAS/IML Software and Matrix Computations

Factor analysis using maximum likelihood estimation

stuti — Sat, 10 Dec 2016 07:26:05 GMT

I want to know the program running behind factor analysis maximum likelihood extraction method. Can i get the inbuilt algorithm for that so that i can verify my result mannually?

my program is given bellow:

start mlfa;

s = r_samp;/* r_samp is my generated correlation matrix data*/

/* objective function module */

start objfun (x) global(s, k, p, xpsy, sigma, lambda, omega, f);

xpsy = diag(1/sqrt(x));

A = xpsy*s*xpsy;

Omega = (eigvec(A)) [,1:k];

theta = (eigval(A)) [1:k];

lambda = diag(sqrt(x))*omega*diag(sqrt(abs(theta-1)));

sigma = lambda*t(lambda)+diag(x);

f = trace(inv(sigma)*s)- log((inv(sigma)*s)- p;

return(f);

finish;

/* Gradient function module */

start gradfun(x) global(s, sigma);

invpsy = inv(diag(x));

g = (vecdiag(invpsy*(sigma-s)*invpsy))`;

return(g);

finish;

/* starting values for iteration */

p = ncol(s);

x = log(((1-k/(2#p))#(1/vecdiag(inv(s))));

/* Set up Options and Constraints */

option = {0 0 . 4};

con = j(1,p, 0.00001)//j(1,p,.);

call nlpqn(rc, xr, "objfun", x, option, con) grd = "gradfun";

r_star = lambda*lambda`;

h_vec = vecdiag(r_star);

all_e_vec = eigvec(R_star);

e_vec = all_e_vec[,1:k];

finish;

start mlfa_extract;

run mlfa;

if rc < 0 then m_rc0 = 1; else m_rc0 = 0;

if rc = 3 then m_rc3 = 1; else m_rc3 = 0;

if rc = 6 then m_rc6 = 1; else m_rc6 = 0;

quit;

Is it correct or i need to do any change for mle?

Re: Factor analysis using maximum likelihood estimation

Rick_SAS — Sat, 10 Dec 2016 11:43:32 GMT

A few comments:

1. The discussion forum usually doesn't say whether a program is correct. If you are having a programming problem, we help you with that problem. Is there a particular part of the program that you want help with?

2. Do you have a reference for the computations you are attempting? It is hard to assess a program without knowing the mathematical formulas it is trying to solve.

3. Do you have sample data that we can use?

4. I would eliminate the MLFA module and just put the data definitions and the call to NLPQN as statements in the main program.

5. The MLFA_EXTRACT module seems to be unfinished.

Overall it looks like you are off to a good start. If you haven't read the article "Maximum likelihood estimation in SAS/IML," that is a good article to read.

Re: Factor analysis using maximum likelihood estimation

stuti — Mon, 12 Dec 2016 06:33:02 GMT

Sir,

I am doing ML estimation of factor analysis. I am attaching ibm spss calculation for ML in factor analysis.

I am running my program on manipulated data having 10 variables for samplesize 30 and pre assumed existance of 2 factors. and if i compare my results with sas university edition inbuilt task program of factor analysis then i am getting different results.

so i need to correct my mle function and for that i need help.

Re: Factor analysis using maximum likelihood estimation

stuti — Mon, 12 Dec 2016 06:36:43 GMT

I am also attaching excel file as the sample data.

On that i ran my program and sas factor analysis using mle.

but the answers are different.

Re: Factor analysis using maximum likelihood estimation

Rick_SAS — Mon, 12 Dec 2016 12:56:05 GMT

The tasks in the SAS University Edition generate SAS code, probably by submitting code for PROC FACTOR.

We need to see the SAS code to know what options you have selected for the algorithm because there are many options in PROC FACTOR, and eash option leads to different results.

Thus the first step should be to look at the SAS code that is generated by the task. Click on the CODE tab to see the code. Then look at the PROC FACTOR documentation to determine how the options affect the factor analysis.

It seems like this isn't really a SAS/IML question. It is really a question about understanding what PROC FACTOR code is being submitted when you use the "Factor Analysis" task in SAS UE.

Re: Factor analysis using maximum likelihood estimation

stuti — Tue, 13 Dec 2016 06:20:18 GMT

Yes sir

I am trying to use the same algorithm which is used by sas UE in proc factor method= ml.

but even in sas documentation, i am unable to find ml method steps.

so if anyone can give me proc factor method=ml inbuilt calculation steps (software inbuilt calculation algorithm) then it would be a great help. because I am developing my simulated data and on that I can't call sas UE proc factor option directly, so i need to go for all sas statements which produces proc factor method= ml output.

Re: Factor analysis using maximum likelihood estimation

Rick_SAS — Tue, 13 Dec 2016 11:24:36 GMT

You say

> even in sas documentation, i am unable to find ml method steps

The SAS documentation says that method=ML "performs maximum likelihood factor analysis with an algorithm due to Fuller (1987)" and gives the reference Fuller, W. A. (1987). Measurement Error Models. New York: John Wiley & Sons.

You say

> I am developing my simulated data and on that I can't call sas UE proc factor option directly

Yes, you can. It is easy to simulate data and call PROC FACTOR. Look for my many blog articles on simulating data or my book Simulating Data with SAS (2013) for many examples of calling SAS procedres on simulated data. Good ones to start with are

Re: Factor analysis using maximum likelihood estimation

Rick_SAS — Tue, 13 Dec 2016 15:18:18 GMT

> if anyone can give me proc factor method=ml inbuilt calculation steps (software inbuilt calculation algorithm) then it would be a great help.

Just to be clear, there is no such thing. To run the algorithm you run PROC FACTOR. Since you are posting this question to the SAS./IML group, I will point out that you can call SAS procedures from SAS/IML programs.

Re: Factor analysis using maximum likelihood estimation

stuti — Sun, 18 Dec 2016 07:18:05 GMT

Thank you so much sir.

I am able to use now proc factor in my simulation program.

One more question can i give name as say matrix 'samp' to only the outcome of mle extraction - rotated factor pattern matrix?

and how can i use and save that samp in further program?

once again thanks.

Re: Factor analysis using maximum likelihood estimation

Rick_SAS — Sun, 18 Dec 2016 12:33:45 GMT

I'm sorry, but I do not understand what you are asking. Can you provide a simple example?

In SAS/IML, you can specify the names of SAS data sets at run time. You use optional arguments in the SUBMIT statement to pass data set names to PROC FACTOR.

Re: Factor analysis using maximum likelihood estimation

stuti — Mon, 19 Dec 2016 03:24:03 GMT

Hello sir

I am attachong one output file which contains one highlited rotated pattern matrix.

In the program after calculating till this step i nedd to calculate Bias which is A1(the first matrix in that file) minus that rotated matrix and then i have to go for further calculation in same manner on rotated pattern matrix.

So I just need rotated pattern matrix('Samp') in further calculation out of all other ml outcomes. the ml procedure which I have used in the program is as follows.

proc factor data=WORK.data1 method=ml rotate=varimax noint nfactors=2
norm=kaiser heywood plots=(scree) outstat=WORK.Factor_ml1;
var col1 col2 col3 col4 col5 col6 col7 col8 col9 col10;
run;

How to just save the rotteted pattern like A1. so that i can find bias. Even I am attaching one excel sheet how to find mannually bias and its average.

Re: Factor analysis using maximum likelihood estimation

Rick_SAS — Mon, 19 Dec 2016 13:25:15 GMT

You can use the ODS OUTPUT statement to create a SAS data set from any output from any SAS procedure. The steps are

1) Put ODS TRACE ON; before you run PROC FACTOR.

2) Run PROC FACTOR. In the SAS lof you will see the names of all ODS tables that were created. The one that you want appears to be OrthRotFactPat.

3) Put ODS OUTPUT OrthRotFactPat=RotFactPattern; before the PROC FACTOR step. Rerun. The data set RotFactPattern contains the information you want. You can use PROC PRINT and PROC CONTENTS to figure out names of variables.

Altogether, it will look like this example:

ods trace on;
proc factor data=sashelp.cars method=ml rotate=varimax noint nfactors=2
                 norm=kaiser heywood plots=(scree) outstat=WORK.Factor_ml1;
/* use all numeric vars */
run;
ods trace off;

/* from the SAS log, it looks like the name of the table 
   we want is OrthRotFactPat */
ods output OrthRotFactPat=RotFactPattern;
proc factor data=sashelp.cars method=ml rotate=varimax noint nfactors=2
                 norm=kaiser heywood plots=(scree) outstat=WORK.Factor_ml1;
run;

proc print data=RotFactPattern; run;

Re: Factor analysis using maximum likelihood estimation

stuti — Tue, 20 Dec 2016 02:53:28 GMT

Thank you so much sir

I am able to use Ods trace in my program and getting that rotated factor matrix in my output.

further can you please guide me how to find Bias which is difference between A1 and rotated factor matrix.

In the last post I have attached mannual excel file for this. after calculating rotated factor pattern matix for the generated sample, I need to compare the original population factor matrix A1 and Mle method generated factor matix using comparision parameters like Bias, Mean square Error, Phi, Factor loading sensitivity and specificity. so if I am able to do calculation for bias may be I will try for other parameters in same manner.

Re: Factor analysis using maximum likelihood estimation

Rick_SAS — Tue, 20 Dec 2016 11:00:13 GMT

If the population parameters are in the matrix A1 in a SSA/IML program, then

1) Read the RotFactPattern matrix into the same SAS/IML session.

use RotFactPattern;

read all var {Factor1 Factor2} into R; / * matrix R contains rotated factors */

close RotFactPattern;

2) Compute the statistics. For example, if bias is the different between the parameters and the estimates, then

Bias = A1 - R;

I feel confident that if you work at this that you will succeed. Best wishes and good luck to you.

Re: Factor analysis using maximum likelihood estimation

stuti — Sat, 24 Dec 2016 07:24:38 GMT

Hello sir

I have first developed the population matrix. Then I derived sample correlation matrix and using CREATE converted into work.samp.data and used proc factor process by using submit and endsubmit statements on that data.. Only rotated factor matrix is used as work.samp_ml.data and on that I used different comparison parameters of rotated matrix and pop matrix A1.

Now I want to replicate the same process say for 1000 sample correlation matrices and then I want to take up average of all comparision parameters.

While doing the replication in proc iml I am getting following error:

ERROR: You cannot open WORK.SAMP_ML.DATA for output access with member-level control because WORK.SAMP_ML.DATA is in use by you in resource environment IML (2).

statement : SUBMIT at line....column 1

Can you help me to solve this?

Re: Factor analysis using maximum likelihood estimation

Rick_SAS — Sat, 24 Dec 2016 10:55:09 GMT

That error message will appear if you open a data set, do not close it, and then try to open it again. Use the CLOSE statement, like this:

do i = 1 to 10;

use samp_ml; /* open the data set */

read all var _num_ into X;

CLOSE samp_ml; /* <== this is the statement you need */

end;

Re: Factor analysis using maximum likelihood estimation

stuti — Sun, 25 Dec 2016 06:38:40 GMT

I have already used close statement but in log it gives the message " Cannot close WORK.SMAP_ML; it is not open."

Re: Factor analysis using maximum likelihood estimation

Rick_SAS — Sun, 25 Dec 2016 11:46:14 GMT

If the error is coming from the SUBMIT statement in which you call PROC FACTOR, then you need to close the data set before the SUBMIT statement. PROC FACTOR can't write to a data set that is already in use.

Re: Factor analysis using maximum likelihood estimation

stuti — Sun, 25 Dec 2016 13:00:10 GMT

yes sir may be error was because of submit statement. I have used close statement earlier to submit but for repetition still without any error message in log, I am not getting answers of all calculations.

For your reference I am coping the entire program here, may be little lengthy but if in first line I use rep=1 then I am getting all results but if I change it even to 2 then I don't get the results of Accuracy, bias and all other parameters. So how to do this repetitions?

Thank you so much for your help. I have learned sas programming just because of this and still trying to learn because of the help provided by you.

SAS program:

%let rep=1; * N of samples from each population;

proc iml;

**********************************************************************;

*************** Begin Specification of Data Conditions ***************;

**********************************************************************;

p = 10; *20; *40; *60;

k = 2; *2; *4; *8;

d_frac = 1; *.05; *.25; *.50; *.75; *.95;

Commun_type = 1;

N_pops = 1; * N of populations to generate;

nn1 = 100000;

means=j(1,p,0);

variance = j(1,p,1);

********************************************************************;

start Make_PopR(nvars,nfactors,commun_type,A1,R);

* program to generate the population covariance matrices

* Inputs are nvars = number of variables in the matrix

nfactors = number of factors

commun = Type of communality

Output is R = population correlation matrix;

bp1=j(nvars,nvars,0);

bp2=uniform(bp1);

if commun_type=3 then do;

bp3=(bp2*.2999999)+.55;

b1square=round(diag(bp3),.1);

end;

if commun_type=1 then do;

bp3=(bp2*.2999999)+.15;

b1square=round(diag(bp3),.1);

end;

if commun_type=2 then do;

bp3=(bp2*.6999999)+.15;

b1square=round(diag(bp3),.1);

end;

B1=b1square##.5;

b3square=I(nvars)-b1square;

B3=b3square##.5;

A1tilde1=j(nvars,nfactors,0);

A1tilde2=round((uniform(A1tilde1)*(nfactors-.00000001))-.5);

A1tilde=A1tilde2;

do j=2 to nfactors;

do i=1 to nvars;

if j<nfactors then do;

A1tilde[i,j]=round(((nfactors-.00000001-sum(A1tilde[i,1:j-1]))

*uniform(0))-.5);

end;

if j=nfactors then do;

A1tilde[i,nfactors]=nfactors-sum(A1tilde[i,1:nfactors-1])-1;

end;

x=normal(A1tilde);

x2=x##2;

d=j(nvars,nfactors,0);

do j=1 to nfactors;

do i=1 to nvars;

d[i,j]=(sum(x2[i,1:nfactors]))##-.5;

end;

cvec=j(1,nfactors,0);

do j=1 to nfactors;

cvec[1,j]=round((uniform(0)*.2999999)+.65,.1);

end;

c=j(nvars,1,1)*cvec;

c2=c##2;

ones=j(nvars,nfactors,1);

y=A1tilde#c + d#x#((ones-c2)##.5);

k=.2;

z=j(nvars,nfactors,0);

do j=1 to nfactors;

do i=1 to nvars;

z[i,j]=((1+k)*y[i,j]*(y[i,j]+abs(y[i,j])+k))/((2+k)*(abs(y[i,j])+k));

end;

z2=z##2;

g=j(nvars,nfactors,0);

do j=1 to nfactors;

do i=1 to nvars;

g[i,j]=(sum(z2[i,1:nfactors]))##-.5;

end;

A1star=g#z;

A1=B1*A1star;

A3star=I(nvars);

A3=B3*A3star;

R=A1*A1`+A3*A3`;

Finish;

start gendata2a(NN1,seed1,variance,bb,cc,dd,mu,r_matrix,YY,p,d_frac);

L = eigval(r_matrix);

neg_eigval = 0;

do r = 1 to nrow(L);

if L[r,1] < 0 then neg_eigval = 1;

end;

if neg_eigval = 0 then do; * matrix is positive definite, so use the Cholesky root approach;

COLS = NCOL(r_matrix);

G = ROOT(r_matrix);

YY=rannor(repeat(seed1,nn1,COLS));

YY = YY*G;

do r = 1 to NN1;

do c = 1 to COLS;

YY[r,c] = (-1*cc) + (bb*YY[r,c]) + (cc*YY[r,c]##2) + (dd*YY[r,c]##3);

YY[r,c] = (YY[r,c] * SQRT(variance[1,c])) + mu[1,c];

end;

if neg_eigval = 1 then do; * matrix is not positive definite, so use the PCA approach;

COLS = NCOL(r_matrix);

V = eigvec(r_matrix);

do i = 1 to nrow(L);

do j = 1 to ncol(V);

if L[i,1] > 0 then V[j,i] = V[j,i] # sqrt(L[i,1]);

if L[i,1] <= 0 then V[j,i] = V[j,i] # sqrt(.000000001);

end;

YY=rannor(repeat(seed1,nn1,COLS));

YY = V*YY`;

YY = YY`;

do r = 1 to NN1;

do c = 1 to COLS;

YY[r,c] = (-1*cc) + (bb*YY[r,c]) + (cc*YY[r,c]##2) + (dd*YY[r,c]##3);

YY[r,c] = (YY[r,c] * SQRT(variance[1,c])) + mu[1,c];

end;

do r = 1 to nn1;

do c = 1 to (p*d_frac);

if yy[r,c] < 0 then yy[r,c] = 0;

else if yy[r,c] = 0 then yy[r,c] = 1;

else if yy[r,c] > 0 then yy[r,c] = 1;

end;

finish;

start gendata2b(NN2,seed1,variance,bb,cc,dd,mu,r_matrix,YY,p,d_frac);

L = eigval(r_matrix);

neg_eigval = 0;

do r = 1 to nrow(L);

if L[r,1] < 0 then neg_eigval = 1;

end;

if neg_eigval = 0 then do; * matrix is positive definite, so use the Cholesky root approach;

COLS = NCOL(r_matrix);

G = ROOT(r_matrix);

YY=rannor(repeat(seed1,nn2,COLS));

YY = YY*G;

do r = 1 to NN2;

do c = 1 to COLS;

YY[r,c] = (-1*cc) + (bb*YY[r,c]) + (cc*YY[r,c]##2) + (dd*YY[r,c]##3);

YY[r,c] = (YY[r,c] * SQRT(variance[1,c])) + mu[1,c];

end;

if neg_eigval = 1 then do;

COLS = NCOL(r_matrix);

V = eigvec(r_matrix);

do i = 1 to nrow(L);

do j = 1 to ncol(V);

if L[i,1] > 0 then V[j,i] = V[j,i] # sqrt(L[i,1]);

if L[i,1] <= 0 then V[j,i] = V[j,i] # sqrt(.000000001);

end;

YY=rannor(repeat(seed1,nn2,COLS));

YY = V*YY`;

YY = YY`;

do r = 1 to NN2;

do c = 1 to COLS;

YY[r,c] = (-1*cc) + (bb*YY[r,c]) + (cc*YY[r,c]##2) + (dd*YY[r,c]##3);

YY[r,c] = (YY[r,c] * SQRT(variance[1,c])) + mu[1,c];

end;

do r = 1 to nn2;

do c = 1 to (p*d_frac);

if yy[r,c] < 0 then yy[r,c] = 0;

else if yy[r,c] = 0 then yy[r,c] = 1;

else if yy[r,c] > 0 then yy[r,c] = 1;

end;

finish;

* +----------------------------------------------------------------+

Subroutine computes bias, MSE and accuracy of factor loadings

Inputs: L = matrix of sample pattern coefficients

Lambda = matrix of population pattern coefficients

* +-----------------------------------------------------------------+;

START ACCURACY(L,Lambda,OK_Load,OK_NoLoad,bias_loadings,MSE_loadings,

pattern_accuracy,perfect_accuracy,

pattern_accuracy30,perfect_accuracy30);

p = nrow(L);

k = ncol(L);

pop_loaded = J(p,1,0);

do variable = 1 to p;

do factor = 1 to k;

if abs(Lambda[variable,factor]) >= .30 then pop_loaded[variable,1] = 1;

end;

sample_loaded = J(p,1,0);

do variable = 1 to p;

do factor = 1 to k;

if abs(L[variable,factor]) >= .30 then sample_loaded[variable,1] = 1;

end;

OK_Load = 0;

OK_NoLoad = 0;

do variable = 1 to p;

if sample_loaded[variable,1] = 1 & pop_loaded[variable,1] = 1 then

OK_Load = OK_Load + 1;

if sample_loaded[variable,1] = 0 & pop_loaded[variable,1] = 0 then

OK_NoLoad = OK_NoLoad + 1;

end;

OK_Load = OK_Load / sum(pop_loaded); * proportion of variables with loading agreement;

if sum(pop_loaded) = p then OK_NoLoad = .;

if sum(pop_loaded) < p then OK_NoLoad = OK_NoLoad/(p-sum(pop_loaded));

* +------------------------------------------------+

Compute bias and MS Error in the factor loadings

* +------------------------------------------------+;

bias_loadings=j(p,k,0);

MSE_loadings=j(p,k,0);

do row =1 to p;

bias_loadings[row,]=L[row,]-lambda[row,];

MSE_loadings[row,]=(L[row,]-lambda[row,])##2;

end;

* +--------------------------------------------+

Compute pattern accuracy using the .30 thumb

* +--------------------------------------------+;

pattern_accuracy = j(p,k,0);

perfect_accuracy = 0;

pattern_accuracy30 = j(p,1,0);

perfect_accuracy30 = 0;

do variable = 1 to p;

NoLoad = 0;

do factor = 1 to k;

if (abs(Lambda[variable,factor]) >= .30 &

abs(L[variable,factor]) >= .30)

then do;

pattern_accuracy[variable,factor] = 1;

pattern_accuracy30[variable,1] = 1;

end;

if (abs(Lambda[variable,factor]) < .30 &

abs(L[variable,factor]) < .30) then do;

pattern_accuracy[variable,factor] = 1;

NoLoad = NoLoad + 1;

end;

if NoLoad = k then pattern_accuracy30[variable,1] = 1;

end;

if mean(pattern_accuracy) = 1 then perfect_accuracy = 1;

if mean(pattern_accuracy30) = 1 then perfect_accuracy30 = 1;

finish;

*+_________________________________________________________________+

Subroutine to calculate coefficient of congruence - mean phi

+_________________________________________________________________+;

start get_phi(lambda,L,meanphi);

k = ncol(lambda);

num = vecdiag(lambda`*L);

den1 = vecdiag(lambda`*lambda);

den2 = vecdiag(L`*L);

phi_k = J(k,1,0);

meanphi = 0;

n_terms = 0;

do m = 1 to k;

if (den1[m,1]>0 & den2[m,1]>0) then do;

phi_k[m,1] = num[m,1] / sqrt(den1[m,1]#den2[m,1]);

meanphi = meanphi + phi_k[m,1];

n_terms = n_terms + 1;

end;

if (den1[m,1] = 0 | den2[m,1] = 0) then do;

phi_k[m,1] = 999;

end;

meanphi = meanphi / n_terms;

finish;

start factor_scores(Lambda,L,sampdat,R_fscores,Non_zero);

*+________________________________________________________________________+

Compute Factor Score Estimates

+_______________________________________________________________________+;

p = nrow(L);

k = ncol(L);

SC =J(p,k,0);

SCP =J(p,k,0);

Non_zero = J(k,1,0);

do j = 1 to p;

do kk = 1 to k;

if L [j,kk] >= .3 then SC[j,kk]=1;

if L [j,kk] <= -.3 then SC[j,kk]=-1;

if Lambda [j,kk] >= .3 then do;

SCP[j,kk]=1;

Non_zero[kk,1] = 1;

end;

if Lambda [j,kk] <= -.3 then do;

SCP[j,kk]=-1;

Non_zero[kk,1] = 1;

end;

FSE = sampdat*SC;

FS = sampdat*SCP;

FS_FSE = FSE||FS;

R_matrix = corr(FS_FSE);

R_matrix2 = R_matrix[k+1:2#k,1:k];

R_fscores = vecdiag(R_matrix2);

finish;

Do pop_num = 1 to N_pops; * Loop for 10 populations;

run Make_PopR(p,k,commun_type,A1,R_pop);

Lambda = A1;

numr = r_pop[+,+] - p;

deno = r_pop[+,+];

ratio = numr/deno;

f2_pop = (p/(p-1))*ratio;

r2_pop = f2_pop/(1+f2_pop);

corr = r_pop;

seed1=round(1000000*rannor(0));

*nn = 100000;

chg = 1;

cycle = 0;

corr_tmp = corr;

do until (chg = 0);

run gendata2a(NN1,seed1,variance,1,0,0,means,corr_tmp,sim_data,p,d_frac);

sim_corr = corr(sim_data);

resid_m = sim_corr - corr;

tot_res = sum(abs(resid_m));

if cycle = 0 then do;

best_corr = corr_tmp;

best_res = tot_res;

end;

if cycle > 0 then do;

if tot_res < best_res then do;

best_corr = corr_tmp;

best_res = tot_res;

end;

if tot_res < (.005#(((p-1)#p)/2)) then CHG = 0; * Convergence!;

if cycle > 30 then do;

if tot_res < (.01#(((p-1)#p)/2)) then CHG = 0; * Convergence!;

end;

if cycle > 200 then CHG = 0;

if CHG = 1 then corr_tmp = corr_tmp - resid_m; * adjust template and simulate another

large sample;

cycle = cycle + 1;

if CHG = 0 then do;

end;

do S_Size =1 to 1; * Loop for sample sizes;

if S_Size = 1 then Sampsize2=3*P;

if S_Size = 2 then Sampsize2=2*p;

if S_Size = 3 then Sampsize2=p;

if S_Size = 4 then Sampsize2=4*p;

if S_Size=5 then Sampsize2=5*p;

* Loop for 1000 Samples;

**********************************************************:

* Begin data generation *;

**********************************************************;

seed1=round(1000000*ranuni(0));

nn2 = sampsize2;

corr_tmp = best_corr;

run gendata2b(NN2,seed1,variance,1,0,0,means,corr_tmp,sim_data,p,d_frac);

sampdat = sim_data;

r_samp=corr(sampdat);

print sampdat;

print A1;

CREATE SAMP FROM SAMPDAT ;APPEND FROM SAMPDAT; CLOSE SAMP;

close work.smap_ml;

SUBMIT;

ods output OrthRotFactPat=samp_ml;

proc factor data=WORK.samp method=ml priors=smc rotate=varimax noint nfactors=2

norm=kaiser heywood ;

var col1 col2 col3 col4 col5 col6 col7 col8 col9 col10;

run;

proc print data=samp_ml;

run;

ENDSUBMIT;

use work.samp_ml; /* open the data set */

read all var {"Factor1" "Factor2" }into L;

close work.samp_ml;

run ACCURACY(L,Lambda,OK_Load,OK_NoLoad,bias_loadings,MSE_loadings,

pattern_accuracy,perfect_accuracy,

pattern_accuracy30,perfect_accuracy30);

run get_phi(lambda,L,meanphi);

run factor_scores(Lambda,L,sampdat,R_fscores,Non_zero);

if &rep = 1 then do;

_ok_load_m = ok_load;

_ok_Noload_m = ok_Noload;

_Bias_Loadings_m = Bias_Loadings;

_MSE_Loadings_m = MSE_Loadings;

_Pattern_Accuracy_m = Pattern_Accuracy;

_Perfect_Accuracy_m = Perfect_Accuracy;

_Pattern_Accuracy30_m = Pattern_Accuracy30;

_Perfect_Accuracy30_m = Perfect_Accuracy30;

_meanphi_m = meanphi;

_R_Fscores_m = R_Fscores;

_Non_zero_m = Non_zero;

if &rep > 1 then do;

_ok_load_m = _ok_load_m + ok_load;

_ok_Noload_m = _ok_Noload_m + ok_Noload;

_Bias_Loadings_m = _Bias_Loadings_m + Bias_Loadings;

_MSE_Loadings_m = _MSE_Loadings_m + MSE_Loadings;

_Pattern_Accuracy_m = _Pattern_Accuracy_m + Pattern_Accuracy;

_Perfect_Accuracy_m = _Perfect_Accuracy_m + Perfect_Accuracy;

_Pattern_Accuracy30_m = _Pattern_Accuracy30_m + Pattern_Accuracy30;

_Perfect_Accuracy30_m = _Perfect_Accuracy30_m + Perfect_Accuracy30;

_meanphi_m = _meanphi_m + meanphi;

_R_Fscores_m = _R_Fscores_m + R_Fscores;

_Non_zero_m = _Non_zero_m + Non_zero;

end;

if &rep = 1 then nsamples = 1;

if &rep > 1 then nsamples = nsamples +1;

*End Replications Loop*;

_ok_load_m = _ok_load_m/NSamples;

_ok_Noload_m = _ok_Noload_m/NSamples;

_Bias_Loadings_m = _BIas_Loadings_m/NSamples;

_MSE_Loadings_m = _MSE_Loadings_m/NSamples;

_Pattern_Accuracy_m = _Pattern_Accuracy_m/NSamples;

_Perfect_Accuracy_m = _Perfect_Accuracy_m/NSamples;

_Pattern_Accuracy30_m = _Pattern_Accuracy30_m/NSamples;

_Perfect_Accuracy30_m = _Perfect_Accuracy30_m/NSamples;

_meanphi_m = _meanphi_m/NSamples;

_R_Fscores_m = _R_Fscores_m/NSamples;

_Non_zero_m = _Non_zero_m/NSamples;

print Sampsize2 commun_type d_frac k p;

print _OK_Load_m _OK_NoLoad_m

_Perfect_Accuracy30_m _Perfect_Accuracy_m _Meanphi_m;

print _Pattern_Accuracy30_m _R_Fscores_m _Non_zero_m;

print _Bias_Loadings_m;

print _MSE_Loadings_m;

print _Pattern_Accuracy_m;

end; end; end;

*End Sample Size Loop*;

*End loop for 10 populations*;

quit;

Re: Factor analysis using maximum likelihood estimation

stuti — Sun, 25 Dec 2016 14:58:14 GMT

Able to modify my program.. now getting results for even repetation...

Thank you so much sir for all the guidence given to me...

still I have to do some modifications in program but may be now I am able to do it.

May be repetation is time consuming but still will be trying to improve it.

If in future would have any problem will come back once again

thank you sir.