NB. multivar.ijs for Multivariate Analysis of Dispersion
NB. Solves Y = XB where Y = n by k, X = n by p and B= p by k.
NB. Design matrices with unequal sample size factors are solved by subtraction. 
NB. Independent effects can be solved directly.
NB. So far:
NB. mulfac implementation- 
NB. (a)- allows entry of Y data matrix and X design matrix. 
NB. (b)- prints out the X'X matrix with columns numbered 0-basis
NB. (c)- prompts for how many effects are to be measured directly
NB. (d)- prompts for the name and columns of each effect.
NB. (e)- prompts for how many effects are to be measured by subtraction.  
NB. (f)- prompts for the name and columns of each effect.
NB. (g)- prints out the ADISP table with Chi Sq and F tests with df.
NB. (h)- deals with special cases of P = 1 or 2 AND Q = 1 or 2
NB. (i)- outputs a Boxed matrix of the MARG matrix and the factor covariances.
NB. (j)- MARG matrix contains the factor names and df.
NB. (k)- add J version printing
NB. (l)- add compression of upper triangular version of the SS matrix
NB. (m)- add formatting of the statistics and SS matrix
NB. (n)- revert to old redu because solution by subtraction did not work
NB.
NB. adinfo implementation-
NB. (a)- inputs boxed output of mulfac.
NB. (b)- extracts MARG & factorial and residual covariance matrices.
NB. (c)- extracts the dimensions & isolates the residual covariance matrix.
NB. (d)- does a test of additional info.
NB.
NB. useful implementation-
NB. (a)- loads in and parses boxed output of mulfac
NB. (b)- problems with matrix dimensions when a depth of 1 is involved.
NB. (c)- Control structures in useful lines[18,31] needed changes for J4.05d involving 'end.'
NB. (d)- 'useful' is not working in J4.06a
NB. 
NB.   For a test, run this script.  It will generate two test Y matrices
NB.     and two test X matrices:  Yab, Yabab, Xab and Xabab.
NB.   To do a test-run, start by invoking the program:

NB.   Rtest =: mulfac 'Test run.'
NB.
NB.   Out =: adinfo Rtest
NB.
NB.    ... this will elicit queries for data input of a Y matrix and an 
NB.   X matrix.  Respond with Yab and Xab, or Yabab and Xab respectively.
NB.
NB.   Mulfac will respond by solving the factorial matrix equation Y = Xb
NB.    for b and then query one for the number of factors to be tested 
NB.    directly and the names of factors and their design column indices.
NB.   After the last factor index is entered, mulfac will ask how many
NB.    factors are to be tested by subtraction.  As above the factor names 
NB.    and indices are queried.
NB.   An ANOVA table is created after the last factor index is entered.
NB.   Then adinfo can be enlisted to take Rtest and process further.

NB. Here are the current J routines:

NB. Default title
title=: 'Default'

NB. How about: QQ=: ([: (1!:1) 1:) [ ] (1!:2) 2:
NB. answer=. 22}. ask 'enter something here: '
ask=: 3 : 0
NB. append line to ijx window (keep trailing blanks)
wd'smappend "',y.,'" 0'
1!:1[1 NB. read from keyboard (ijx window)
)

NB. Find the last non blank in y.
last_blank=: # - 1: #.~ ] = ' '"_

NB. Matrix operations mp, trans, inv
load'c:\j405\user\projects\jgk\matrix_primitives.ijs'

NB. Steele and Torrie test data Table 11-2
load'c:\j405\user\projects\jgk\sandt-tab11-2.ijs'

NB. Steele and Torrie test data Table 15-12 (loads Ygui and Xgui)
load'c:\j405\user\projects\jgk\sandt-tab15-12.ijs'

NB. generation of a random column 2 of data from Steele & Torrie data
Yabab=:Yab,.Yab+ (20 1$30?30)-15

Yababab=: Yabab,. Yabab + 0.25*(20 2$40?40)

NB. Solution of the general linear model
lreg=: 4 : 0
yty=: (trans y.) mp y.
xtx=: (trans x.) mp x.
xty=: (trans x.) mp y.
B=: (inv xtx)mp xty
)


NB. Main factorial regression program, mulfac.
mulfac=: 3 :0
  Title=: y.
 MESS=:' '
 9!:17 (1) NB. Set box columns to centered.
 ] 9!:14''
 Y =: ". 16}. ask 'Enter Y matrix: '
 X =: ". 16}. ask 'Enter X matrix: '
 RDF=: (,1 {.$Y) - (,1}.$X)
 junk =. X lreg Y
 junk=. (3 3$('Solve Y = XB'; ' ';'X''X';'B';'Index' ;( 3j0 ": trans z); B;(z=:(IK,1)$i.IK=:1{.$xtx);xtx)) 1!:2 (2)
 junk=. ask ''
 R0 =: yty - REG=: RED i.IK
 R=: (1,$yty)$0
 MARG=:1 20$' '
 L=: 1 3 3 $0
 ne =: ". 30}. ask 'Number to be tested directly? '
 ne SOU 'RV=: RED V2'
 ne =: ". 36}. ask 'Number to be tested by subtraction? '
 J=: 1{.$xtx
 ne SOU 'RV=: REG -(trans B)mp XR mp B=:(Vc getxy xty)inv XR=:(Vc=. V2 redu i. J) getxx xtx '
 junk=. (Title,' on ',9{.": 6!:0 '' ) 1!:2 (2)
 junk=. ('Analysis of Dispersion Table printed at: ',9}.": 6!:0 '' ) 1!:2 (2)
 MARG=:1 0}.MARG
 L=: 1 0 0}.L
 junk=. (2 2$ 'Source          df';'Lambda         ChiSq  df   n   p   q df1  df2  F-ratio'; MARG;(10j8 10j3 4j0 4j0 4j0 4j0 4j0 5j1 9j3 ":((1{.$L),9)$,L)) 1!:2 (2)
 R=: 1 0 0}.R,R0
 junk=. (2 2$ 'Source          df';(,11j0 3j0 ":>clabels);(MARG=: MARG, 1 20$'Residual SSy.x  ', 4j0 ":RDF);(14j3 ": uindex utri R)) 1!:2 (2)
 junk=. (2 2$ 'Source          df';'Ry.x';MARG;((1{.$R),(*/1}.$R))$,R)
)

RED=: 3 :0
V=: y.
red=.(trans V{xty)mp V{B
)

NB. Computation of Lambda and appending of Chi-Sq and F stats with FLAM
LAMB=: 4 :0
R0=. x.
RV=. y.
lambda=: (det R0)% det R0 + RV
N=: RDF + Q=: $,V2
P=: 1}.$R0
if. (lambda > 0) do. goto_L2. else. goto_L1. end.
label_L1. MESS=: MESS, 'lambda negative, '
lambda=: 1
label_L2. lambda=: 3 3 $(_1 }. lambda,(-M*^.lambda),(P*Q),N,P,Q,M=: N - 0.5 * P+Q+1), FLAM lambda
)

NB. List of credits and services to be offered
LOGO=: 3 :0
a=. 'Multivariate Extension of the Generalised Linear Model: Y=XB'
b=. 'written by Joe Kunkel in APL (1980), translated to J (1999)'
c=. 'r =: mulfac ''title'' = solution of Y=XB for B and Ry.x'
d=. 't =: adinfo r = multivariate test for aditional information'
e=. 'u =: useful r = test of usefullness of multivariate covariates'
f=. 'based on C.R.Rao (1965) Linear Statistical Inference and its Applications'
junk=. (6 1 $ a;b;c;d;e;f) 1!:2 (2)
)

NB. List of credits and services to be offered
TRY=: 3 :0
a=. 'Multivariate Extension of the Generalised Linear Model: Y=XB'
b=. 'Steele and Torrie test data Table 11-2'
c=. 'Yabab=:Yab,.Yab  and  Yababab=: Yabab,. Yabab + 0.25*(20 2$40?40)'
d=. 'Steele and Torrie test data Table 15-12 (loads Ygui and Xgui)'
e=. 'Ygui;Xgui:';Ygui;Xgui
f=. 'Try using them as input once you have invoked Rout=: mulfac ''title'''
g=. ('kernel32 GetComputerNameA i *c *i' 15!:0 ((30#' ');,30))
h=. 'J-version: ',(9!:14'')
junk=. (8 1 $ a;b;c;d;e;f;g;h) 1!:2 (2)
)


NB. query and apply calculation of independent or subtraction effects
SOU=: 4 : 0
NE=. x.
FUNC=. y.
if. (NE > 0) do. goto_S0. else. goto_S2. end.
label_S0. I=: 0
label_S1. I=: I + 1
SOUR =: 9}. ask 'Source ',(":I),' '
SOURCE =. 1 15{.(1,rs=.$SOUR)$ SOUR
V2=: ,". (rs+19)}. ask 'Column indices of ',SOUR,'? '
MARG=: MARG, SOURCE,. 1 5$ 5j0":$V2 
". FUNC
L=: L, (R0 LAMB RV)
R=: R, RV
if. (NE > I) do. goto_S1. else. goto_S2. end.
label_S2. return.
)

NB. calculation of F plus its df for general P,Q and special cases.
FLAM=: 3 : 0
L=. ,y.
t=. N
NB. Separate cases where Q = 1 or 2.
if. (Q < 3) do. goto_F1. else. goto_F3. end.
label_F1. if. (Q=2) do. goto_F2. else. F=:(1-L)*(t-P)%(L*P)
D1=: P
D2=: t-P
goto_F8. end.
label_F2. F=:(1-L^0.5)*(t-P)%(P*L^0.5)
D1=: 2*P
D2=: 2*(t-P-1)
goto_F8.
NB. Separate cases where P = 1 or 2.
label_F3. if. (P<3) do. goto_F4. else. goto_F6. end.
label_F4. if. (P=2) do. goto_F5. else. F=:(1-L)*(t-Q)%(L*Q)
D1=: Q
D2=: t-Q
goto_F8. end.
label_F5. F=:(1-L^0.5)*(t-Q)%(Q*L^0.5)
D1=: 2*Q
D2=: 2*(t-Q-1)
goto_F8. 

label_F6. M=: N - 0.5*P+Q+1
S=: ((_4+D1*D1=:P*Q)% (_5+(P*P)+Q*Q))^0.5 NB. made M, D1 and S global 7/20/99
LAM=. 0.25 * D1-2
F=: ,((1-L)*D2=: (M*S)-2*LAM)%D1*L=.L^(%S) NB. made D2 global 7/20/99
label_F8. F=: D1,D2,F
)

NB. redu=: ([:-.e.)#[ NB. reduce index y. by elements of x.

redu=: 4 : 0  NB. reduce index y. by elements of x.
ind=. y.
knot=. $ not=. x.
if. (knot < 2) do. (not=. 1$not) end.
not=. 1 - +/not =/ ind
ind=. _1+ not * ind +1
ind =. /:~ ind
drop=. +/ _1 = ind
ind=. drop }. ind
)


getxx=: 4 : 0  NB. reduce matrix y. using index x. 
matrix =. y.
index =. x.
matrix =. index { matrix
matrix =. trans index { trans matrix
)

getxy=: 4 : 0  NB. reduce matrix y. using index x. 
matrix =. y.
index =. x.
matrix =. index { matrix
)

NB. Test for Additional information { adinfo boxed(Marg|Residual) }
adinfo=: 3 :0
k=:1{.j=.2{. 2}.$MARG=:> _1 1{. y.
width=: 1}.j
MARG=: (k,width)$,MARG
SOUR=: 0 _5 }.MARG
label_L0. 
P=: (_1 {. $R=: > _1 _1 {. y.)^0.5
R=:(k,P,P)$, R
R0=: (P,P)$, ((k-1),0 0)}. R
a=.'Indexed Sources:'
junk=. (1 16 $ a) 1!:2 (2)
9!:17 (1) NB. Set box columns to centered.
junk=. (2 2$ 'Index';'Source            df';(2j0 ":(k,1)$i.k);MARG) 1!:2 (2)
RDF=: ". _1 _5 {. MARG
junk=. ask ''
I =: ". 44}. ask 'In what source are you interested?<0 exits> '
if. (I=0) do. goto_L99. else.
OUTPUT=: 1 48$' '
lb=. 1 + last_blank source=: , _1 15{.((I+1),15){.SOUR
source=: lb {. source
RI=:(P,P)$ , (1,P,P){. (I,0 0)}. R
junk=. (('total df= ', 4j0 ": T=: RDF+Q),' and factor ',source,' df= ', 4j0 ": Q=: ". _1 _5 {.(((1-k-I),0)}.MARG)) 1!:2 (2)
V2=: Q $ 1
junk=. ask ''
label_L1. 
RJ=: $ J=:  ". A=: 47}. ask 'Given trait indices:<998=new source, 999=print> '
if. (J=998) do. goto_L0. elseif. (J=999) do. goto_L10. elseif. do.
Ps=:P
R0b=: J getxx R0
RIb=: J getxx RI
LU=: 1 1{. R0b LAMB RIb
junk=. ('Lambda= ',(, 8j6 ": LU),'   of interest: ',3j0 ": J) 1!:2 (2)
junk=. ask ''
label_L2. 
PS=: $ K=:  ". B=: 57}. ask 'Is there additional info in the traits:<999=new base set> '
if. (K=999) do. goto_L1. else.
OUT=:30 $' '
R0k=:(Vk=. K redu J) getxx R0
RIk=: Vk getxx RI
L=: LU% 1 1{. R0k LAMB RIk
N=: T - RJ - P
OUT=: 24 $ B,'/',(":Vk),OUT
F=: FLAM L
OUTPUT=:OUTPUT, 1 48$OUT,(,6j0 ": 2{.F),(,12j4 ": 2}.F)
goto_L2.
end.
end.
end.
label_L10.
junk=. ('')1!:2 (2)
 a=. (Title,': for Factor: ',source )
 b=. (": 6!:0)
 b=. ('Test of Additional Information')
 c=. ('Traits/Covariates            DFu     DFl      F ')
 d=. (1 0 }. OUTPUT)
junk=. (4 1 $ a;b;c;d) 1!:2 (2)
goto_L0.
label_L99.
LOGO 1
) 


NB. Test for usefulness of correlation corrections { useful boxed(Marg|Residual) }
useful=: 3 :0
k=:1{.j=.2{. 2}.$MARG=:> _1 1{. y.
width=: 1}.j
MARG=: (k,width)$,MARG
SOUR=: 0 _5 }.MARG
P=: (_1 {. $R=: > _1 _1 {. y.)^0.5
R=:(k,P,P)$, R
R=: (P,P)$, ((k-1),0 0)}. R
Rpad=:(R,0),.0
corequire 'jwatch'
w2=: 'Rpad' conew 'jwatch'
a=.'Indexed Sources:'
junk=. (1 16 $ a) 1!:2 (2)
9!:17 (1) NB. Set box columns to centered.
junk=. ask ' '
junk=. (2 2$ 'Index';'Source            df';(2j0 ":(k,1)$i.k);MARG) 1!:2 (2)
RDF0=: ". _1 _5 {. MARG
junk=. ask ' '
V2=: ". A=: 37}. ask 'Give covariate variables:<999=print> '
if. ((1{.V2)=999) do. goto_L1. elseif. (COVdf=:(1{.$V2)=0) do. goto_L11.  elseif. do. goto_L12. end.
label_L11. COVdf=:1
label_L12. RDF=: RDF0 - COVdf
junk=. ask ' '
J=: ". B=: 26}. ask 'Give dependent variables: '
R0s=: $R0=: J getxx R
junk=. (1 $$R0s ) 1!:2 (2)
if. ((1{.$J)=0) do. goto_L13.  else. goto_L14. end.
label_L13. R0=: 1 1 $R0
label_L14.
R0pad=:(R0,0),.0
corequire 'jwatch'
w2=: 'R0pad' conew 'jwatch'
if. 1=1 do. junk=.1 end.
R11=: V2 getxx R
if. ((1{.$V2)=0) do. goto_L15.  else. goto_L16. end.
label_L15. R11=: 1 1 $R11
label_L16.
R11pad=:(R11,0),.0
corequire 'jwatch'
w2=: 'R11pad' conew 'jwatch'

R12=: trans V2 getxy (J getxy R)   NB. R12=: trans V2 getxy trans (J getxy R)
if. ((1{.$R12)=0) do. goto_L17. elseif.((1{.$$R12)=1) do. goto_L18. elseif. do. goto_L19. end.
label_L17. R12=: 1 1 $R12
label_L18. nv=:,1 {. $R12
R12=: (1,nv)$R12
label_L19.
R12pad=:(R12,0),.0
corequire 'jwatch'
w2=: 'R12pad' conew 'jwatch'
shape =:$RP=: R12 mp ((inv R11)mp  trans R12) NB. R0=:R0-RP=:  R12 mp ((inv R11)mp trans R12)
R0=:(shape $ R0)-RP
if. ((1{.$J)=0) do. goto_L20.  else. goto_L21. end.
label_L20. R0=: 1 1 $R0
label_L21.
R0pad2=:(R0,0),.0
corequire 'jwatch'
w2=: 'R0pad2' conew 'jwatch'

F0 =: _1 3 {. R0 LAMB RP

label_L1.
)

NB. getting an upper triangular representation

utri =: 3 : 0
sh =. i. $ y.
sh =. (sh = 1|:sh)+ (sh < 1|:sh)
sh * y.
)

NB. Cartesian product
  CP=: {@(,&<)

NB. Harvest the rectangular form of the base upper triangular SS matrices
uindex =: 3 : 0
indy=.1+ 1 {. sh =. i. $ y.
nm=. 1{.$y.
ns=. _1{.$y.
sh =. ,indy* 1{.(sh = 1|:sh)+ (sh < 1|:sh =. sh +1)
indy=. (nm,ns*ns)$,y.
clabels=: ,(i. ns) CP i.ns
clabels=: (+/i.ns)}."1 clabels /:"1 sh
(+/i.ns)}."1 indy /:"1 sh
)


TRY ''

NB.    Rtest =: mulfac 'Test run.'
NB. Enter Y matrix: Yab
NB. Enter X matrix: Xab
NB. ÚÄÄÄÄÄÄÂÄÄÄÄÄÂÄÄÄÄÄÄÄÄÄÄÄ¿
NB. ³      ³     ³    X'X    ³
NB. ÃÄÄÄÄÄÄÅÄÄÄÄÄÅÄÄÄÄÄÄÄÄÄÄÄ´
NB. ³  B   ³Index³  0 1 2 3  ³
NB. ÃÄÄÄÄÄÄÅÄÄÄÄÄÅÄÄÄÄÄÄÄÄÄÄÄ´
NB. ³24.246³  0  ³20  0  0  0³
NB. ³ 7.927³  1  ³ 0 20  0  0³
NB. ³ 3.041³  2  ³ 0  0 10  0³
NB. ³_4.361³  3  ³ 0  0  0 10³
NB. ÀÄÄÄÄÄÄÁÄÄÄÄÄÁÄÄÄÄÄÄÄÄÄÄÄÙ

NB. Number to be tested directly? 4
NB. Source 1 mean
NB. Column indices of mean? 0
NB. Source 2 a
NB. Column indices of a? 1
NB. Source 3 b
NB. Column indices of b? 2
NB. Source 4 ab
NB. Column indices of ab? 3
NB. Number to be tested by subtraction? 0
NB. Test run. on 2000 9 28
NB. Analysis of Dispersion Table printed at:  15 59 5.123
NB. ÚÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÂÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄ¿
NB. ³ Source          df ³ Lambda ChiSq df n p q df1 df2 F-ratio ³
NB. ÃÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÅÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄ´
NB. ³mean               1³0.0313021 53.6931 1 17 1 1 1 16 495.147³
NB. ³a                  1³ 0.232132  22.637 1 17 1 1 1 16 52.9263³
NB. ³b                  1³  0.80424 3.37678 1 17 1 1 1 16 3.89455³
NB. ³ab                 1³ 0.666408 6.29073 1 17 1 1 1 16 8.00933³
NB. ÀÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÁÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÙ
NB. ÚÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÂÄÄÄÄÄÄÄ¿
NB. ³ Source          df ³ Ry.x  ³
NB. ÃÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÅÄÄÄÄÄÄÄ´
NB. ³mean               1³11757.4³
NB. ³a                  1³1256.75³
NB. ³b                  1³92.4768³
NB. ³ab                 1³190.183³
NB. ³Residual SSy.x    16³379.923³
NB. ÀÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÁÄÄÄÄÄÄÄÙ