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)- changed file address unix path for Mac
NB. (l)- changed ask function to remove windows dependent call.
NB. (m)- added plotting of B matrix.
NB. (n)- added XLABEL deffinitions.
NB. (o)- remove mean as a plotted effect.
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. 
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. answer=. 22}. ask 'enter something here: '
ask=: 3 : 0
NB. append line to ijx window (keep trailing blanks)
NB. wd'smappend "',y.,'" 0'
junk=. (y.)1!:2 (2)
1!:1[1 NB. read from keyboard (ijx window)
)

load 'plot'

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

NB. Matrix operations mp, trans, inv
load'\Users\jgkunkel\j501b\user\projects\jgk\matrix_primitives.ijs'

NB. Steele and Torrie test data Table 11-2
load'\Users\jgkunkel\j501b\user\projects\jgk\sandt-tab11-2.ijs'

NB. Steele and Torrie test data Table 15-12 (loads Ygui and Xgui)
load'\Users\jgkunkel\j501b\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.
 Y =: ". 0}. ask 'Enter Y matrix: '
 X =: ". 0}. ask 'Enter X matrix: '
 RDF=: (,1 {.$Y) - (,1}.$X)
 junk =. X lreg Y
 junk=. (3 3$(' '; ' ';'X''X';'B';'Index' ;(trans z); B;(z=:(IK,1)$i.IK=:1{.$xtx);xtx)) 1!:2 (2)
 junk=. ask ''
 R0 =: yty - REG=: RED i.IK
 BSAVE=:B
 R=: (1,$yty)$0
 MARG=:1 20$' '
 L=: 1 3 3 $0
 ne =: ". 0}. ask 'Number to be tested directly? '
 ne SOU 'RV=: RED V2'
 xlab1=: xlab
 ne =: ". 0}. 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 '
 xlab2=: xlab
 junk=. (Title,' on ',10{.": 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;(((1{.$L),9)$,L)) 1!:2 (2)
 R=: 1 0 0}.R,R0
 junk=. (2 2$ 'Source          df';'Ry.x';(MARG=: MARG, 1 20$'Residual SSy.x  ', 4j0 ":RDF);((1{.$R),(*/1}.$R))$,R) 1!:2 (2)
 do_plot''
 
)

do_plot=: 3 : 0
NB. initialize form here
pd 'reset'
redb=. 1 0}.BSAVE
nplots=. , _1{. $redb
max=. >./,redb
min=. <./,redb
span=. <.(max-min)%10
tic=. <.span%4

NB. dat =. 1}., (0 ,-(nplots-1)) }. BSAVE
dat =. trans redb
pd 'type bar'
pd 'ytic ', ": span, tic
pd 'xlabel ',xlab1
pd 'title ',Title
pd dat
pd 'show'
)



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=. ,'Operating System: ', 1{.(10 -(dum=.9!:12''))}. SYS
junk=. (7 1 $ a;b;c;d;e;f;g) 1!:2 (2)
)

 SYS=: 9 9 $ 'PC       PC-386   Windows  MacintoshOS2      Unix     Win95/98/Win CE   Other     '



NB. query and apply calculation of independent or subtraction effects
SOU=: 4 : 0
xlab=:''
NE=. x.
FUNC=. y.
if. (NE > 0) do. goto_S0. else. goto_S2. end.
label_S0. I=: 0
label_S1. I=: I + 1
SOUR =: 0}. ask 'Source ',(":I),' '
SOURCE =. 1 15{.(1,rs=.$SOUR)$ SOUR
V2=: ,". 0}. ask 'Column indices of ',SOUR,'? '
if. (I = 1) do. goto_L1. else. goto_L2. end.
label_L2. xlab=: xlab, make_quoted V2
label_L1. 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.
)

 
make_quoted=: 3 : 0
ndex=: , $y.
nstr=: , $ SOUR
qbl=: (ndex, 1)$'"'
qstr=: (ndex,nstr)$ SOUR
qndex=: (ndex, 1)$y.
Xlabels =: ,qbl,.qstr,. (":qndex),. (ndex,2)$'" '
)


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
)

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

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 =: ". 0}. 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=: 0}. 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=: 0}. 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=: 0}. 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=: 0}. 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.
)

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. ??????????????????????????????