writing formula

[dustychair]

Hello,

I wrote a formula in SAS IML. In that formula there are two type of parameters. One of them is from a existing file and the other one is created by SAS. When I run the program it works but the results are weird. I want to be sure that my code is correct. Could you check my code please. The formula which I work is attached.

data parms;
	INFILE "C:\mirt\parms.txt";
		INPUT a1 a2 a3 d1;

run;

proc iml;

  call randseed(0);

  use parms;

  read  all var {a1 a2 a3} into a;

  read all var {d1} into d;

   mean = {0,0,0};  *change these values to transfrom theta1 and theta2 to a different mean;

   cov = {1 0.45 0.45,

            0.45 1 0.45,

0.45 0.45 1};  *change the covariance to modify the correlation between dimensions;

  theta = RANDNORMAL(5000,mean,cov);

T=theta`;

z=j(5000,30);

p=j(5000,30);

u=j(5000,30);

  do j=1 to 5000;

        do i=1 to 30;

			p[j,i]= exp((a[i,]*T[,j])+d[i,])/(1+(exp((a[i,]*T[,j])+d[i,])));
				u[j,i]=ranuni(0);

                if p[j,i]<u[j,i] then z[j,i]=0;

                else if p[j,i]>=u[j,i]then z[j,i]=1;

        end;

end;

create response from z;

append from z;

names1={theta1 theta2 theta3};

create theta from theta [colname=names1];

append from theta;

quit;

data keep2;

  set response;

  file "C:\mirt\mirt_data.txt";

  put (col1-col30)(1.0);

run;

proc corr data=keep2 alpha;

  var col1-col30;

  ods output cronbachalpha=alpha;

  run;

data alpha_need;

  set alpha;

  if _n_=1;

  keep alpha;

  file "c:\mirt\alpha\alpha.txt" mod;

  put alpha;

run;

data theta;

  set theta;

  file "C:\mirt\theta.txt";

  put theta1 theta2 theta3;

run;

[Rick_SAS]

I don't know what "the results are weird" means.  What are you seeing that is different than you expect?

Your loop appears to reconstruct the formula in the paper.  I don't see anything obviously wrong, other than the code being slightly inefficient.

I suggest that you debug and check the code on a smaller example, like 5x3 instead of 5000x30.  That will enable you to check the computations by hand to make sure they are correct.

[Rick_SAS]

The program is much easier to understand (and more efficient) if you vectorize the computations instead of using two DO loops.  In vector form, the main computation requires two calls:

p = logistic( theta*A` + d` );

call randgen(z,"Bernoulli", p);

The computation of P is equivalent to the longer computation that you have:

do j=1 to 5000;

   do i=1 to 30;

      p[j,i]= exp((a[i,]*T[,j])+d[i,])/(1+(exp((a[i,]*T[,j])+d[i,])));

   end;

end;

[dustychair]

Dr. Wicklin,

Thank you SO much for your answer. I really appreciate it. I will work on both method. By the way "weird results" means different results from literature. Now, I'm pretty sure these results are  because of my design.