topic Re: writing formula in SAS/IML Software and Matrix Computations

writing formula

dustychair — Thu, 05 Jun 2014 16:14:16 GMT

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;

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;

Re: writing formula

Rick_SAS — Fri, 06 Jun 2014 15:43:12 GMT

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.

Re: writing formula

Rick_SAS — Mon, 09 Jun 2014 11:22:05 GMT

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;

Re: writing formula

dustychair — Mon, 09 Jun 2014 17:27:16 GMT

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.