***  NOTE this is a subroutine.


* *************************** NOTE   NOTE *********************************.
*  1) No more than 10 independent variables.
*  2) Standardize the independents ( good thing to do in any regression analysis).
*  3) Look at options under edit and set work space to 9999.
 * 4) This may take several minutes to run.
 * 5) If this runs, then run allpos_2  - this make the file manageable and produces some nice graphs.
* 6) you data file WILL BE REPLACED. So save it.



*MACRO CALL.
*allpos Y x1 x2   - xP.

allpos  mpg engine weight horse accel.

* Now Run the following.

allpos_2.



DEFINE allpos(!positional !TOKENS(1) /!positional !cmdend).
 SET MXLOOP =9999.
MATRIX.
GET W/VARIABLES = !1 !2 /Names=varnames/FILE=* /Missing=Omit.

COMPUTE NC = NCOL(W).
compute nam = make (1,nc,"a").
loop i1 = 2 to nc .
compute nam(i1) = varnames(i1).
end loop.
compute nam(1) ="const".

 Print nc.
  do if (nc <12) .
compute ni = nc-1.
print ni.
compute ntot=2**ni -1.
* Print ntot.
compute all1= make(ntot,5,0).
compute var1 =make(ntot,ni+1,0).
loop i1 = 1 to ntot.
compute var1(i1,1) = 1.
end loop.
loop i1 = 1 to ni.
compute var1(1,i1+1) =1.
end loop.

 Print ntot.
COMPUTE NR =NROW(W).

loop ll1 = 1 to ntot.
compute all1(ll1,5) = ll1.
end loop.

COMPUTE y=MAKE(NR,1,0).
COMPUTE BIGX =MAKE(NR,NC,1).
LOOP II = 1 TO NR.
COMPUTE y(II,1)=W(II,1).
LOOP JJ  = 1 TO NC-1.
COMPUTE BIGX(II,JJ+1)=W(II,JJ+1).
END LOOP.
END LOOP.
compute jmatrix =make(nr,1,1).
compute meany =t(y)*jmatrix/nr.
compute a1 = y -meany*jmatrix.
compute ssto = t(a1)*a1.
* Print ssto.


* * Print BIGX.
* * Print x.
compute a = bigx*inv(T(bigx)*bigx)*t(bigx).
compute mset = t(y-a*y)*(y-a*y)/(nr-nc).
* Print mset.
compute ssrt = ssto-(nr-nc)*mset.
* Print ssrt.
compute sse = mset*(nr-nc).
* Print sse.
COMPUTE ADJRSQ = 1 - ((NR-1)/(NR - NI -1))*SSE/SSTO.
* Print ADJRSQ.
COMPUTE  Mallows = ni+1.
* Print Mallows.
Compute cont = 1.
compute all1(cont,1) = mallows.
compute all1(cont,2) = adjrsq.
compute all1(cont,3) = mset.
compute all1(cont,4) = ni+1.
do if (ni >1).
* subsets of size 1.
Compute SUBx = make(nr,2,1).

loop i = 1 to ni . 
* print i.

loop k = 1 to nr.
compute subx(k,2) =bigx(k,i+1).

end  loop.
* Print i.
* Print j.
  * Print subx.
compute a = subx*inv(T(subx)*subx)*t(subx).
compute mse = t(y-a*y)*(y-a*y)/(nr-2).
* Print mse1.
compute sse = (nr-2)*mse.
COMPUTE ADJRSQ = 1 - ((NR-1)/(NR - 2))*(SSE/SSTO).
* Print ADJRSQ1.
compute mallows = (sse/mset) -(nr -2*2).
* Print mallows.
Compute cont = cont+1.
* Print cont.
 compute all1(cont,1) = mallows.
compute all1(cont,2) = adjrsq.
compute all1(cont,3) = mse.
compute all1(cont,4) = 2.
compute var1(cont,i+1) = 1.
 * Print cont.
* print i.
end loop.
end if.


do if (ni >2).
* Subsets of Size 2.
Compute SUBx = make(nr,3,1).

loop i = 1 to ni - 1. 
loop j = i+1 to ni.
* print j.
loop k = 1 to nr.
* print k.
compute subx(k,2) =bigx(k,i+1).
compute subx(k,3) = bigx(k,  j+1).
end  loop.
*  Print i.
*  Print j.
* Print subx.
compute a = subx*inv(T(subx)*subx)*t(subx).
compute mse = t(y-a*y)*(y-a*y)/(nr-3).
* Print mse2.
compute sse = (nr-3)*mse.
COMPUTE ADJRSQ = 1 - ((NR-1)/(NR - 3))*(SSE/SSTO).
* Print ADJRSQ2.
compute mallows = (sse/mset) -(nr -2*3).
* Print mallows.
* print ni.
 * print cont.
Compute cont = cont+1.
compute all1(cont,1) = mallows.
compute all1(cont,2) = adjrsq.
compute all1(cont,3) = mse.
compute all1(cont,4) = 3.
compute var1(cont,i+1) = 1.
compute var1(cont,j+1) = 1.
end loop.
end loop.
end if.

do if (ni >3).
* Subsets of Size 3.
Compute SUBx = make(nr,4,1).

loop i = 1 to ni - 2. 
loop j = i+1 to ni-1.
loop j2 = j+1 to ni.

loop k = 1 to nr.
compute subx(k,2) =bigx(k,i+1).
compute subx(k,3) = bigx(k,j+1).
compute subx(k,4) = bigx(k,j2+1).
end  loop.
* print subx.
* Print i.
*  Print j.
* Print subx.
compute a = subx*inv(T(subx)*subx)*t(subx).
compute mse = t(y-a*y)*(y-a*y)/(nr-4).
* Print mse2.
compute sse = (nr-4)*mse.
COMPUTE ADJRSQ = 1 - ((NR-1)/(NR - 4))*(SSE/SSTO).
* Print ADJRSQ2.
compute mallows = (sse/mset) -(nr -2*4).
* Print mallows.

Compute cont = cont+1.
compute all1(cont,1) = mallows.
compute all1(cont,2) = adjrsq.
compute all1(cont,3) = mse.
compute all1(cont,4) = 4.
compute var1(cont,i+1) = 1.
compute var1(cont,j+1) = 1.
compute var1(cont,j2+1) = 1.
end loop.
end loop.
end loop.
end if.
print ni.
do if (ni>4).

* subsets of size 4.

Compute SUBx = make(nr,5,1).

loop i = 1 to ni - 3. 
loop j = i+1 to ni-2.
loop j2 = j+1 to ni-1.
loop j3 = j2+1 to ni.

loop k = 1 to nr.
compute subx(k,2) =bigx(k,i+1).
compute subx(k,3) = bigx(k,j+1).
compute subx(k,4) = bigx(k,j2+1).
compute subx(k,5) = bigx(k,j3+1).
end  loop.
* print subx.
* Print i.
*  Print j.
* Print subx.
compute a = subx*inv(T(subx)*subx)*t(subx).
compute mse = t(y-a*y)*(y-a*y)/(nr-5).
* Print mse2.
compute sse = (nr-5)*mse.
COMPUTE ADJRSQ = 1 - ((NR-1)/(NR - 5))*(SSE/SSTO).
* Print ADJRSQ2.
compute mallows = (sse/mset) -(nr -2*5).
* Print mallows.

Compute cont = cont+1.
compute all1(cont,1) = mallows.
compute all1(cont,2) = adjrsq.
compute all1(cont,3) = mse.
compute all1(cont,4) = 5.
compute var1(cont,i+1) = 1.
compute var1(cont,j+1) = 1.
compute var1(cont,j2+1) = 1.
compute var1(cont,j3+1) = 1.

end loop.
end loop.
end loop.
end loop.
end if.


* subset of size 5.
do if (ni >5).

Compute SUBx = make(nr,6,1).

loop i = 1 to ni - 4. 
loop j = i+1 to ni-3.
loop j2 = j+1 to ni-2.
loop j3 = j2+1 to ni-1.
loop j4 = j3+1 to ni.
loop k = 1 to nr.
compute subx(k,2) =bigx(k,i+1).
compute subx(k,3) = bigx(k,j+1).
compute subx(k,4) = bigx(k,j2+1).
compute subx(k,5) = bigx(k,j3+1).
compute subx(k,6) =bigx(k,j4+1).
end  loop.
* print subx.
* Print i.
*  Print j.
* Print subx.
compute a = subx*inv(T(subx)*subx)*t(subx).
compute mse = t(y-a*y)*(y-a*y)/(nr-6).
* Print mse2.
compute sse = (nr-6)*mse.
COMPUTE ADJRSQ = 1 - ((NR-1)/(NR -6))*(SSE/SSTO).
* Print ADJRSQ2.
compute mallows = (sse/mset) -(nr -2*6).
* Print mallows.

Compute cont = cont+1.
compute all1(cont,1) = mallows.
compute all1(cont,2) = adjrsq.
compute all1(cont,3) = mse.
compute all1(cont,4) = 6.
compute var1(cont,i+1) = 1.
compute var1(cont,j+1) = 1.
compute var1(cont,j2+1) = 1.
compute var1(cont,j3+1) = 1.
compute var1(cont,j4+1) =1.

end loop.
end loop.
end loop.
end loop.
end loop.
end if.



* subset of size 6.
do if (ni >6).

Compute SUBx = make(nr,7,1).

loop i = 1 to ni - 5. 
loop j = i+1 to ni-4.
loop j2 = j+1 to ni-3.
loop j3 = j2+1 to ni-2.
loop j4 = j3+1 to ni-1.
loop j5 = j4+1 to ni.
loop k = 1 to nr.

compute subx(k,2) =bigx(k,i+1).
compute subx(k,3) = bigx(k,j+1).
compute subx(k,4) = bigx(k,j2+1).
compute subx(k,5) = bigx(k,j3+1).
compute subx(k,6) =bigx(k,j4+1).
compute subx(k,7) =bigx(k,j5+1).
end  loop.
* print subx.
* Print i.
*  Print j.
* Print subx.
compute a = subx*inv(T(subx)*subx)*t(subx).
compute mse = t(y-a*y)*(y-a*y)/(nr-7).
* Print mse2.
compute sse = (nr-7)*mse.
COMPUTE ADJRSQ = 1 - ((NR-1)/(NR -7))*(SSE/SSTO).
* Print ADJRSQ2.
compute mallows = (sse/mset) -(nr -2*7).
* Print mallows.

Compute cont = cont+1.
compute all1(cont,1) = mallows.
compute all1(cont,2) = adjrsq.
compute all1(cont,3) = mse.
compute all1(cont,4) = 7.
compute var1(cont,i+1) = 1.
compute var1(cont,j+1) = 1.
compute var1(cont,j2+1) = 1.
compute var1(cont,j3+1) = 1.
compute var1(cont,j4+1) =1.
compute var1(cont,j5+1) =1.

end loop.
end loop.
end loop.
end loop.
end loop.
end loop.
end if.

* subset of size 7.
do if (ni >7).

Compute SUBx = make(nr,8,1).

loop i = 1 to ni - 6. 
loop j = i+1 to ni-5.
loop j2 = j+1 to ni-4.
loop j3 = j2+1 to ni-3.
loop j4 = j3+1 to ni-2.
loop j5 = j4+1 to ni-1.
loop j6 = j5+1 to ni.
loop k = 1 to nr.

compute subx(k,2) =bigx(k,i+1).
compute subx(k,3) = bigx(k,j+1).
compute subx(k,4) = bigx(k,j2+1).
compute subx(k,5) = bigx(k,j3+1).
compute subx(k,6) =bigx(k,j4+1).
compute subx(k,7) =bigx(k,j5+1).
compute subx(k,8)=bigx(k,j6+1).
end  loop.
* print subx.
* Print i.
*  Print j.
* Print subx.
compute a = subx*inv(T(subx)*subx)*t(subx).
compute mse = t(y-a*y)*(y-a*y)/(nr-8).
* Print mse2.
compute sse = (nr-8)*mse.
COMPUTE ADJRSQ = 1 - ((NR-1)/(NR -8))*(SSE/SSTO).
* Print ADJRSQ2.
compute mallows = (sse/mset) -(nr -2*8).
* Print mallows.

Compute cont = cont+1.
compute all1(cont,1) = mallows.
compute all1(cont,2) = adjrsq.
compute all1(cont,3) = mse.
compute all1(cont,4) = 8.
compute var1(cont,i+1) = 1.
compute var1(cont,j+1) = 1.
compute var1(cont,j2+1) = 1.
compute var1(cont,j3+1) = 1.
compute var1(cont,j4+1) =1.
compute var1(cont,j5+1) =1.
compute var1(cont,j6+1) =1.

end loop.
end loop.
end loop.
end loop.
end loop.
end loop.
end loop.
end if.
* subset of size 8.
do if (ni >8).

Compute SUBx = make(nr,9,1).

loop i = 1 to ni - 7. 
loop j = i+1 to ni-6.
loop j2 = j+1 to ni-5.
loop j3 = j2+1 to ni-4.
loop j4 = j3+1 to ni-3.
loop j5 = j4+1 to ni-2.
loop j6 = j5+1 to ni-1.
loop j7 = j6+1 to ni.
loop k = 1 to nr.

compute subx(k,2) =bigx(k,i+1).
compute subx(k,3) = bigx(k,j+1).
compute subx(k,4) = bigx(k,j2+1).
compute subx(k,5) = bigx(k,j3+1).
compute subx(k,6) =bigx(k,j4+1).
compute subx(k,7) =bigx(k,j5+1).
compute subx(k,8)=bigx(k,j6+1).
compute subx(k,9)=bigx(k,j7+1).
end  loop.
* print subx.
* Print i.
*  Print j.
* Print subx.
compute a = subx*inv(T(subx)*subx)*t(subx).
compute mse = t(y-a*y)*(y-a*y)/(nr-9).
* Print mse2.
compute sse = (nr-9)*mse.
COMPUTE ADJRSQ = 1 - ((NR-1)/(NR -9))*(SSE/SSTO).
* Print ADJRSQ2.
compute mallows = (sse/mset) -(nr -2*9).
* Print mallows.

Compute cont = cont+1.
compute all1(cont,1) = mallows.
compute all1(cont,2) = adjrsq.
compute all1(cont,3) = mse.
compute all1(cont,4) = 9.
compute var1(cont,i+1) = 1.
compute var1(cont,j+1) = 1.
compute var1(cont,j2+1) = 1.
compute var1(cont,j3+1) = 1.
compute var1(cont,j4+1) =1.
compute var1(cont,j5+1) =1.
compute var1(cont,j6+1) =1.
compute var1(cont,j7+1) =1.


end loop.
end loop.
end loop.
end loop.
end loop.
end loop.
end loop.
end loop.
end if.

* subset of size 9.
do if (ni >9).

Compute SUBx = make(nr,10,1).

loop i = 1 to ni - 8. 
loop j = i+1 to ni-7.
loop j2 = j+1 to ni-6.
loop j3 = j2+1 to ni-5.
loop j4 = j3+1 to ni-4.
loop j5 = j4+1 to ni-3.
loop j6 = j5+1 to ni-2.
loop j7 = j6+1 to ni-1.
loop j8 = j7+1 to ni.
loop k = 1 to nr.

compute subx(k,2) =bigx(k,i+1).
compute subx(k,3) = bigx(k,j+1).
compute subx(k,4) = bigx(k,j2+1).
compute subx(k,5) = bigx(k,j3+1).
compute subx(k,6) =bigx(k,j4+1).
compute subx(k,7) =bigx(k,j5+1).
compute subx(k,8)=bigx(k,j6+1).
compute subx(k,9)=bigx(k,j7+1).
compute subx(k,10) = bigx(k,j8+1).
end  loop.
* print subx.
* Print i.
*  Print j.
* Print subx.
compute a = subx*inv(T(subx)*subx)*t(subx).
compute mse = t(y-a*y)*(y-a*y)/(nr-10).
* Print mse2.
compute sse = (nr-10)*mse.
COMPUTE ADJRSQ = 1 - ((NR-1)/(NR -10))*(SSE/SSTO).
* Print ADJRSQ2.
compute mallows = (sse/mset) -(nr -2*10).
* Print mallows.

Compute cont = cont+1.
compute all1(cont,1) = mallows.
compute all1(cont,2) = adjrsq.
compute all1(cont,3) = mse.
compute all1(cont,4) = 10.
compute var1(cont,i+1) = 1.
compute var1(cont,j+1) = 1.
compute var1(cont,j2+1) = 1.
compute var1(cont,j3+1) = 1.
compute var1(cont,j4+1) =1.
compute var1(cont,j5+1) =1.
compute var1(cont,j6+1) =1.
compute var1(cont,j7+1) =1.
compute var1(cont,j8+1) =1.


end loop.
end loop.
end loop.
end loop.
end loop.
end loop.
end loop.
end loop.
end loop.
end if.

 *  print all1.
*print var1.
Compute all2 ={all1,var1}.
compute cnames ={"Mallows","ADJRSQ","MSE","P", "Model"}.
compute can1 ={cnames,nam}.
*print all2.
  save all2 /outfile = */ names =can1. 
 else.
 Print /Title " More Than 10  Independent variables".
 end if.
end matrix.
!ENDDEFINE. 



define allpos_2 ( ).

compute p1 =p.
execute.


SORT CASES BY
  p (a) adjrsq (d) .

compute filter1 = 0.
execute.
Matrix.
Get x/variables = all/file = */names = names.
* print x.
compute nr = nrow(x).
compute nc = ncol(x).
compute x(1,nc) = 1.
compute x(2,nc) = 1.
compute x(nr,nc)=1.
loop i = 3 to nr-1.
do if (x(i,nc-1) <>x(i-1,nc-1)).
compute x(i,nc) = 1.
compute x(i+1,nc) = 1.

end if.
end loop.
save x /outfile = */names = names.
end matrix.


FILTER OFF.
USE ALL.
SELECT IF(filter1=1).
EXECUTE .




GRAPH
  /SCATTERPLOT(OVERLAY)=p p  WITH p1 mallows (PAIR) BY model (IDENTIFY)
  /MISSING=LISTWISE
  /TITLE= "Mallow's Cp Versus P".

GRAPH
  /SCATTERPLOT(BIVAR)=p WITH adjrsq BY model
  /MISSING=LISTWISE
  /TITLE= 'Adjusted Rsq Versus P'.
GRAPH
  /SCATTERPLOT(BIVAR)=p WITH mse BY model
  /MISSING=LISTWISE
  /TITLE= 'MSE Versus P'.


!enddefine.