Re: MATRIX LIMITATIONS, Negative Variance Components etc...

[David Marso]

PLEASE IGNORE THIS PREVIOUS POSTING.
I HAVE NOTED A FEW ISSUES TO BE CORRECTED SOON.
DELETED FROM NABBLE (did not propagate to UGA).
Please accept my apology.
David

=====================
To manage your subscription to SPSSX-L, send a message to
[hidden email] (not to SPSSX-L), with no body text except the
command. To leave the list, send the command
SIGNOFF SPSSX-L
For a list of commands to manage subscriptions, send the command
INFO REFCARD

[David Marso]

* I posted this early yesterday to be alerted by Ryan and Bruce that all was not well in Houston
* or  somewhere (Oz?) ).  DELETED from Nabble, but original still alive at UGA (how to purge that?).
* Issues had to do with :
* 1. DATASETS (Who knows?  Code worked for me, consistently failed for R & B).
* 2. I had Transposed the CHOL and had arrived at incorrect correlations *.
* 3. Pending: OMS issues -see below-

* So, let's try this again *.
* Related threads:
* 1 http://spssx-discussion.1045642.n5.nabble.com/Interpretation-of-a-valid-negative-variance-component-ICC-td5718961.html 
* 2 http://spssx-discussion.1045642.n5.nabble.com/SPSS-MATRIX-LIMITATION-td5719047.html
------------------------------------------------------------------------------------------------------------------
* Inspired by two threads initiated by Ryan Black, I conjured the following.
* Initially I merely desired to generalize the #aij elements to something other than 3x3 (see thread 1).
* Pondering on this I ran a simple matrix program against a 3x3 with r=-.35 with a simple CHOLESKY.
* Sure enough the #aij are elements from a Cholesky decomposition. When I attempted to run this
* against a 4x4 with r=-.35 it failed due to an NPD matrix. Figured probably good to notify user so
* grabbed the eigenvalues, checked the last one for >0. If not EXIT otherwise go for it (notice how that
* works). Moving right along I also decided to add an optional Orthonormalization routine (GACK, would
* prefer to just use MATRIX rather than FACTOR, but just didn't have time to deal with that).

* Finally, my first Foray into OMS (someone please help fix this: I want the MATRIX text to trap
* the NPD * error..., but don't want the notes! I want to totally suppress the FACTOR output but
* couldn't get it to function via OMS so I used my SET RESULTS OFF ERRORS OFF hack to deal with that).

*Re Thread 2.
* 1. Technique for generation of dummy data without INPUT PROGRAM using MATRIX.
* 2. Hand off to Main system to generate X vectors of RV.UNIFORM & Optionally orthonormalize.
* 3. Hand back to MATRIX to do the remaining simulation (implemented using a matrix 'submacro').

Some issues: Things sometimes get FUBAR with the DATASET business. I believe it is now fixed with the NEWFILE DATASET/CLOSE ALL double tap! This could certainly be refined, but I only have so much time! Happy Simulations.

* If you are still awake, here is the code *. ********************************************************************. * Created by David Marso, March 2013 **. * DISCLAIMERS: **. * AS IS! IF THIS DOESN'T WORK FOR YOU THEN FEEL FREE TO FIX IT **. * DO NOT REPOST ANYWHERE WITHOUT WRITTEN PERMISSION (I'm easy) **. ********************************************************************.

**** From Bruce ****.
* David, based on something I learned playing around with other code
* you posted a while ago, I suspected the problem was related to using
* datasets rather than writing to files on the disk.  So I made the
* changes needed to test that idea, and hey presto, it worked.

* After the fact, I thought I should indicate where the edits are,
* so have prefaced them with:  ** BW edits .
* I hope I found them all again.
* THANKS BRUCE YOU ARE A LIFE SAVER!!!!!!!*


*  If you are still awake, here is the code *.
********************************************************************.
* Created by David Marso, March 2013                              **.  
* DISCLAIMERS:                                                    **.
*  AS IS! IF THIS DOESN'T WORK FOR YOU THEN FEEL FREE TO FIX IT   **.
*  DO NOT REPOST ANYWHERE WITHOUT WRITTEN PERMISSION  (I'm easy) **.
********************************************************************.
                     
DEFINE ORTHONORMALIZE (!POS !CHAREND ("/") / !POS !TOKENS(1) ).
** Need Matrix expressions for normalization **.
** In Lieu of that: QAD using FACTOR scores  **.
+  PRESERVE.
+  *******************************************************************.
+  ** Would be nice if I could get OMS to work and lose the Tex-Fix **.
+  *******************************************************************.
+  SET ERRORS OFF RESULTS OFF.
+  FACTOR VAR !1 /CRITERIA FACTORS(!2)/ SAVE REG (ALL,FX).
+  RESTORE.
!ENDDEFINE.

DEFINE GOFORIT (Sigma !TOKENS(1)).
+  ** Transform Orthogonal vectors to correlated **.
+  ** Previously had: T(CHOL(#RMat)) OOPS *.
+  COMPUTE E=!Sigma * X * CHOL(#RMat).
+  COMPUTE Y=MAKE(NROW(X),1,0).

+  LOOP #=1 TO NROW(X).
+    COMPUTE Y(#)=E(#,SUBJECT(#)).
+  END LOOP.
+  ************************************************************************.
+  ** BW edit on next line (swapped out Disk File reference for DATASET) **.
+  ************************************************************************.
+  SAVE {PLOT, SUBJECT, Y} / OUTFILE "SimulationOut.sav"/VARIABLES PLOT SUBJECT Y.
!ENDDEFINE .

DEFINE NegativeVarComponentSimulationUsingMatrix (
           sigma !CHAREND ("/")
         / rho   !CHAREND ("/")
         / nplot !CHAREND ("/")
         / nreps !CHAREND ("/")
         / ORTH  !CHAREND ("/") !DEFAULT ('N')
         / DEBUG !CHAREND ("/") !DEFAULT ('N') ).
+  PRESERVE.

+  SET MXLOOPS 10000000.

************************************************.
** THIS NEEDS HELP                            **.
** Warnings                                   **.
** Command OMS. The subcommand If is repeated **.
** Execution of this command stops            **.
************************************************.

+  !IF (!DEBUG !EQ 'N') !THEN .
+ * OMS.
+  OMS /SELECT ALL /IF COMMANDS=['Matrix'] /IF SUBTYPES =['NOTES']/DESTINATION VIEWER=NO.
+  OMS /SELECT ALL /IF COMMANDS=['Factor'] /DESTINATION VIEWER=NO.
+ * OMS.
+  !IFEND .

+  NEW FILE.
+  DATASET CLOSE ALL.
+  *****************************************************************.
+  ** BW edit on next 7 lines (swapped out Disk File for DATASET) **.
+  *****************************************************************.
+  MATRIX.
+    SAVE (MAKE((!nplot*!nreps),(!nreps+2),0))
        / OUTFILE "Simulation.sav"
        / VARIABLES plot subject X1 TO !CONCAT(X,!nreps).
+  END MATRIX.

********************************************************.
+  GET FILE "Simulation.sav".
+  DATASET NAME Simulation.
+  DATASET ACTIVATE Simulation.
+  COMPUTE subject=MOD($CASENUM,!nreps).
+  COMPUTE plot=TRUNC(($CASENUM-.1)/!NREPS)+1.

+  !LET !VLIST = !CONCAT(' X1 TO X',!nreps) .
+  VECTOR X=!VLIST .

** Generate Normal deviates **.
+  LOOP #=1 TO !nreps.
+    COMPUTE X(#)=RV.NORMAL(0,1).
+  END LOOP.
+  RECODE SUBJECT (0=!nreps).

********************************************************.
+  !IF (!ORTH !EQ 'Y') !THEN.
+    Orthonormalize !VLIST / !nreps.
+    DELETE VARIABLES !VLIST.
+    RENAME VARIABLES (!CONCAT('FX1 TO FX',!nreps,"=", !VLIST)).
+  !IFEND.
********************************************************.

********************************************************.
+  MATRIX.
+    GET X        /FILE * /VARIABLES !VLIST.
+    GET SUBJECT  /FILE * /VARIABLES SUBJECT.
+    GET PLOT     /FILE * /VARIABLES PLOT .

+    COMPUTE #RMat= MAKE(!NReps,!NReps,!rho).
+    CALL SETDIAG(#RMAT,1).

+    !IF (!DEBUG !EQ 'Y' ) !THEN
+      PRINT #RMAT / TITLE "Correlation Matrix" .
+    !IFEND.

+    ***********************************************************.
+    ** Check that submitted CORR is PD prior to calling CHOL **.
+    ***********************************************************.
+    CALL EIGEN(#RMAT,EIGVEC,EIGVAL).
+    DO IF (EIGVAL(NROW(EIGVAL)) GT 0).
+      GoForIt Sigma= !Sigma.
+    ELSE.
+      PRINT / TITLE "Correlation Matrix is NOT Positive Definite" .
+      PRINT T(EIGVAL) / TITLE "Eigenvalues of Correlation Matrix".
+    END IF.

+  ** ADDED TO CONFIRM CORRELATIONS ARE CORRECT **.
+    !IF (!DEBUG !EQ Y ) !THEN
+      SAVE E / OUTFILE "TestEForCorr.sav" / VARIABLES !VLIST .
+    !IFEND .
+  END MATRIX.

+  *****************************************************************.
+  ** BW edit on next 2 lines (swapped out Disk File for DATASET) **.
+  *****************************************************************.
+  GET FILE = "SimulationOut.sav".
+  DATASET NAME SimulationOut.
+  DATASET ACTIVATE SimulationOut.

+  !IF (!DEBUG !EQ 'N') !THEN .
+ * OMS.
+  OMSEND.
+  !IFEND .

+  !IF (!DEBUG !EQ Y ) !THEN
+    GET FILE "TestEForCorr.sav".
+    CORRELATIONS VAR !VList.
+  !IFEND .

+  RESTORE.
!ENDDEFINE.
SET PRINTBACK ON.
NegativeVarComponentSimulationUsingMatrix sigma=1 /rho= -.25 /nplot=50 /nreps 4 /ORTH=Y.

** Append this parameter if you want to see a tour of the sausage factory **.
 /DEBUG=Y.