Turn on suggestions

Auto-suggest helps you quickly narrow down your search results by suggesting possible matches as you type.

Showing results for

Options

- RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Mute
- Printer Friendly Page

🔒 This topic is **solved** and **locked**.
Need further help from the community? Please
sign in and ask a **new** question.

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

Posted 05-16-2020 07:37 PM
(818 views)

Hello, I am working to obtain OLS parameter estimates (with no intercept) for each possible pair of observations within subgroups. Let's say the year 1991 represents a group, then let's say there are two subgroups within 1991; the first subgroup contains 3 observations, which produces 3 possible pairs or combination (1-2, 2-3 and 1-3); and the second subgroup contains 4 observations, which produces 6 possible pairs or combinations (1-2, 1-3, 1-4, 2-3, 2-4 and 3-4). The goal is to generate OLS parameter estimates for a model such as Y = B1 + B2 + e (i.e., two indep vars and no intercept). If this sounds like the Theil-Sen method, you are spot on. However, note that this is the *"multivariate version"*, which requires more than computing the slope for y = mx + b. A fish much larger than I from the academic world has determined that OLS with no intercept and two observations will generate Theil-Sen slopes for B1 and B2 mentioned earlier (a modified r-square is necessary for OLS with no intercept, and will typically be between 99.5% and 100% when using the big fish methodology).. I have the "data step" SAS code from the big fish, but having witnessed the power and efficiency of IML, I set out to improve the lengthy processing time of the "data step" code. The IML code I developed below uses the OLS code from one of Rick Wicklin's blog, along with a loop developed by a SAS Community member. It works at the group level, but I lack the IML knowledge to make it work at the subgroup level. I believe the key is to insert a DO loop within the existing DO loop to capture the dynamics of the subgroups. This will require reconstructing a new X matrix for each possible pair within a subgroup in order to compute the two slope values. If interested or helpful, I can supply the code from the big fish. Thank you for any ideas or suggestions! Rick

```
data untrimmed;
input group lag2cvrank lagcfo_ts cfo_ts lag2cfo_ts lag2acr_ts;
cards;
1991 0 155 175 165 35
1991 0 200 225 250 75
1991 0 75 125 135 65
1991 1 350 375 400 55
1991 1 155 175 165 85
1991 1 200 225 250 100
1991 1 75 125 135 125
;
run;
/* find unique BY-group combinations */
proc freq data=untrimmed;
tables group*lag2cvrank / out=FreqOut;
run;
proc iml;
start regress(XY);
*c = allcomb(nrow(XY),2); /* all "N choose 2" combinations of pairs */
*c_rows=nrow(c);
group = XY[1,4]; /* extract group from XY, to be used as a BY variable later */
lag2cvrank = XY[1,5]; /* extract lag2cvrank from XY, to be used as a BY variable later */
/* Extract x from XY */
X = XY[c[i],{1 2}]; /* extract X from XY */ /* extract pairs, i.e, each combo in C */
/* Extract y from XY */
Y = XY[c[,],3]; /* extract Y from XY */
xpx = x`*x; /*cross-products*/
xpy = x`*y;
/*solve linear system*/
/*Solution 1: compute inverse with INV (inefficient)*/
xpxi = inv(xpx); /*form inverse crossproducts*/
b = xpxi*xpy; /*solve for parameter estimates*/
* Or a better solution ***;
/*Solution 2: compute solution with SOLVE. More efficient*/
b = (solve(xpx, xpy))`; /* solve for parameter estimates*/
t = nrow(XY); /* number of rows in XY */
group_col = group;
lag2cvrank_col = lag2cvrank;
return (b || group_col || lag2cvrank_col);
end;
finish;
/* read the BY groups */
use FreqOut nobs NumGroups;
read all var {group lag2cvrank};
close FreqOut;
use work.untrimmed;
create ts_fcst1_IML var {m m1 group_col lag2cvrank_col};
setin work.untrimmed;
setout ts_fcst1_IML;
inVarNames = {"lag2cfo_ts" "lag2acr_ts" "lagcfo_ts" "group" "lag2cvrank"};
do i = 1 to NumGroups; /* for each BY group */
read all var inVarNames into XY
where(group=(group[i]) & lag2cvrank=(lag2cvrank[i]));
print xy;
/* X contains data for i_th group; analyze it */
c = allcomb(nrow(XY),2); /* all "N choose 2" combinations of pairs */
c_rows=nrow(c);
do i = 1 to c_rows; /* this is my feeble attempt at the new DO loop */
G = regress(XY);
/* extract the columns of the matrix */
m=G[,1]; m1=G[,2]; group_col=G[,3]; lag2cvrank_col=G[,4];
append;
end;
end;
close work.untrimmed;
close ts_fcst1_IML;
quit;
data ts_fcst1_IML;
set ts_fcst1_IML;
proc print;run;
```

1 ACCEPTED SOLUTION

Accepted Solutions

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

I developed the code below to solve this problem. The downside is that IML processing time is infinitely greater than than data-step processing time. Thank you SAS Community! Rick

```
data untrimmed;
input group lag2cvrank lagcfo_ts cfo_ts lag2cfo_ts lag2acr_ts;
cards;
1991 0 155 175 165 35
1991 0 200 225 250 75
1991 0 75 125 135 65
1991 1 350 375 400 55
1991 1 155 175 165 85
1991 1 200 225 250 100
1991 1 75 125 135 125
1992 3 100 250 175 35
1992 3 125 175 250 100
1992 3 200 225 300 175
;
run;
/* find unique BY-group combinations */
proc freq data=untrimmed;
tables group*lag2cvrank / out=FreqOut;
run;
proc iml;
start regress(XY);
*c = allcomb(nrow(XY),2); /* all "N choose 2" combinations of pairs */
*c_rows=nrow(c);
c = allcomb(nrow(XY),2); /* all "N choose 2" combinations of pairs */
c_rows=nrow(c);
b_storage_first={0 0};
do i = 1 to c_rows;
XY_NEW=XY[c[i,],];
group = XY[1,4]; /* extract group from XY, to be used as a BY variable later */
lag2cvrank = XY[1,5]; /* extract lag2cvrank from XY, to be used as a BY variable later */
/* Extract x from XY */
X = XY_NEW[,{1 2}]; /* extract X from XY */ /* extract pairs, i.e, each combo in C */
/* Extract y from XY */
Y = XY_NEW[,3]; /* extract Y from XY */
xpx = x`*x; /*cross-products*/
xpy = x`*y;
/*solve linear system*/
/*Solution 1: compute inverse with INV (inefficient)*/
xpxi = inv(xpx); /*form inverse crossproducts*/
b = xpxi*xpy; /*solve for parameter estimates*/
* Or a better solution ***;
/*Solution 2: compute solution with SOLVE. More efficient*/
b = (solve(xpx, xpy))`; /* solve for parameter estimates*/
t = nrow(XY_NEW); /* number of rows in XY */
group_col = J(c_rows, 1, group);
lag2cvrank_col = J(c_rows, 1, lag2cvrank);
If i=1 then b_storage=b_storage_first+b;
else b_storage=b_storage//b;
print b_storage;
end;
return (b_storage || group_col || lag2cvrank_col);
finish;
/* read the BY groups */
use FreqOut nobs NumGroups;
read all var {group lag2cvrank};
close FreqOut;
use work.untrimmed;
create ts_fcst1_IML var {group_col lag2cvrank_col med_b1 med_b2};*m m1;
setin work.untrimmed;
setout ts_fcst1_IML;
inVarNames = {"lag2cfo_ts" "lag2acr_ts" "lagcfo_ts" "group" "lag2cvrank"};
do i = 1 to NumGroups; /* for each BY group */
read all var inVarNames into XY
where(group=(group[i]) & lag2cvrank=(lag2cvrank[i]));
/* X contains data for i_th group; analyze it */
G = regress(XY);
/* extract the columns of the matrix */
m=G[1,1]; m1=G[1,2]; group_col=G[1,3]; lag2cvrank_col=G[1,4]; /* Save matrix components into variables */
/* Compute median slope values from all pairs for group/lag2cvrank combinations */
med_b1=median(m);
med_b2=median(m1);
append;
end;
close work.untrimmed;
close ts_fcst1_IML;
quit;
data ts_fcst1_IML;
set ts_fcst1_IML;
proc print;run;
```

8 REPLIES 8

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

Since it sounds like this is an academic research project, I encourage you to think about the statistical computation you want to make for each combination from the subgroups. Suppose that cases 1 and 2 are being used for the first subgroup and cases 3 and 4 are being used for the second subgroup. How do you estimate B1 and B2 for those observations? Be sure you completely understand the answer and discuss it with your advisor if you have any questions.

1. Write a function that performs that computation to make sure you really understand it. Make sure it also works when cases 1&2 are used for the first subgroup and cases 1&2 are used for the second subgroup.

2. Now think about how to handle all combinations from the subgroups. Use the ALLCOMB function on each subgroup to generate all pairwise combinations from subgroup 1 and from subgroup 2.

3. Leverage the function you wrote in Step 1 to get all estimates over all pairwise combinations from the subgroups.

4. Try to encapsulate steps 1-3 into a single function. The input to the function is a group. The output is the set of Sen-Theil estimates for the subgroups.

5. After you have successfully completed steps 1-4, you can leverage your existing program to read each group, call the function in Step 4, and write the estimates to a SAS data set.

You are creating a big program, so you must break it into smaller pieces and then implement each small piece as efficiently as possible. In computer science lingo, this is sometimes called top-down design and bottom-up construction. Make sure you understand the underlying math/stats at each step before you begin to write any code.

Good luck!

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

I developed the code below to solve this problem. The downside is that IML processing time is infinitely greater than than data-step processing time. Thank you SAS Community! Rick

```
data untrimmed;
input group lag2cvrank lagcfo_ts cfo_ts lag2cfo_ts lag2acr_ts;
cards;
1991 0 155 175 165 35
1991 0 200 225 250 75
1991 0 75 125 135 65
1991 1 350 375 400 55
1991 1 155 175 165 85
1991 1 200 225 250 100
1991 1 75 125 135 125
1992 3 100 250 175 35
1992 3 125 175 250 100
1992 3 200 225 300 175
;
run;
/* find unique BY-group combinations */
proc freq data=untrimmed;
tables group*lag2cvrank / out=FreqOut;
run;
proc iml;
start regress(XY);
*c = allcomb(nrow(XY),2); /* all "N choose 2" combinations of pairs */
*c_rows=nrow(c);
c = allcomb(nrow(XY),2); /* all "N choose 2" combinations of pairs */
c_rows=nrow(c);
b_storage_first={0 0};
do i = 1 to c_rows;
XY_NEW=XY[c[i,],];
group = XY[1,4]; /* extract group from XY, to be used as a BY variable later */
lag2cvrank = XY[1,5]; /* extract lag2cvrank from XY, to be used as a BY variable later */
/* Extract x from XY */
X = XY_NEW[,{1 2}]; /* extract X from XY */ /* extract pairs, i.e, each combo in C */
/* Extract y from XY */
Y = XY_NEW[,3]; /* extract Y from XY */
xpx = x`*x; /*cross-products*/
xpy = x`*y;
/*solve linear system*/
/*Solution 1: compute inverse with INV (inefficient)*/
xpxi = inv(xpx); /*form inverse crossproducts*/
b = xpxi*xpy; /*solve for parameter estimates*/
* Or a better solution ***;
/*Solution 2: compute solution with SOLVE. More efficient*/
b = (solve(xpx, xpy))`; /* solve for parameter estimates*/
t = nrow(XY_NEW); /* number of rows in XY */
group_col = J(c_rows, 1, group);
lag2cvrank_col = J(c_rows, 1, lag2cvrank);
If i=1 then b_storage=b_storage_first+b;
else b_storage=b_storage//b;
print b_storage;
end;
return (b_storage || group_col || lag2cvrank_col);
finish;
/* read the BY groups */
use FreqOut nobs NumGroups;
read all var {group lag2cvrank};
close FreqOut;
use work.untrimmed;
create ts_fcst1_IML var {group_col lag2cvrank_col med_b1 med_b2};*m m1;
setin work.untrimmed;
setout ts_fcst1_IML;
inVarNames = {"lag2cfo_ts" "lag2acr_ts" "lagcfo_ts" "group" "lag2cvrank"};
do i = 1 to NumGroups; /* for each BY group */
read all var inVarNames into XY
where(group=(group[i]) & lag2cvrank=(lag2cvrank[i]));
/* X contains data for i_th group; analyze it */
G = regress(XY);
/* extract the columns of the matrix */
m=G[1,1]; m1=G[1,2]; group_col=G[1,3]; lag2cvrank_col=G[1,4]; /* Save matrix components into variables */
/* Compute median slope values from all pairs for group/lag2cvrank combinations */
med_b1=median(m);
med_b2=median(m1);
append;
end;
close work.untrimmed;
close ts_fcst1_IML;
quit;
data ts_fcst1_IML;
set ts_fcst1_IML;
proc print;run;
```

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

You say the IML solution is taking a lot longer, but there are probably ways in which it could be made more efficient. It is generally a bad idea to 'grow' a matrix inside a loop with syntax like:

`b_storage=b_storage//b;`

as this makes a new matrix for each iteration, and moves data around in memory unnecessarily. It is better to create b_storage at it final size before the loop, and then write one row at a time.

I am guessing the data step solution you mention uses a short cut to compute regression coefficients directly, since there are only two points and no error. You could do the same in IML to speed things up. For an even faster solution, it could be vectorized to eliminate the loop over the combinations.

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

To learn more about Ian's comment, see the article "Friends don't let friends concatenate results inside a loop."

I'll also mention that you are solving the regression problem twice in every loop. You can delete the unnecessary statements for Solution 1.

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

My suggestion is to declare the storage immediately before the loop like so:

`b_storage = j(c_rows, 2, .);`

That is a matrix with as many rows as combinations, 2 columns and filled with missing values. Once you have the matrix b in the loop, then copy it to the ith row of the storage matrix as follows:

`b_storage[i, ] = b;`

The if/then/else statement can then be deleted, as can the definition of matrix b_storage_first.

Incidentally the statements that define group_col and lag2cvrank_col can be moved outside of the loop, they are being overwritten on every iteration of the loop - this will save a bit more time.

Finally, have you thought what happens if you have two data rows with identical pairs of x-values?

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

**Don't miss out on SAS Innovate - Register now for the FREE Livestream!**

Can't make it to Vegas? No problem! Watch our general sessions LIVE or on-demand starting April 17th. Hear from SAS execs, best-selling author Adam Grant, Hot Ones host Sean Evans, top tech journalist Kara Swisher, AI expert Cassie Kozyrkov, and the mind-blowing dance crew iLuminate! Plus, get access to over 20 breakout sessions.

Multiple Linear Regression in SAS

Learn how to run multiple linear regression models with and without interactions, presented by SAS user Alex Chaplin.

Find more tutorials on the SAS Users YouTube channel.