/* From Cameron & Pauling (1978) Supplemental ascorbate in the supportive treatment of cancer: Reevaluation of prolongation of survival times in terminal human cancer. Procceedings of the National Academy of Science, 75, 4548-4542. Table 1 Alive=0 (no) or 1 (yes) The variables are survival time in days */ data vitC; input type $ case gender $ age y1 x1 y2 x2 alive; d1 = y1-x1; d2 = y2-x2; if type='breast' or type='ovary' then delete; if alive=1 then delete; /* Design Matrix A: n x 7 with overall and 6 */ ao=1; a1=0; a2=0; a3=0; a4=0; a5=0; a6=0; if type='stomach' then a5=1; else if type='bronchus' then a2=1; else if type='colon' then a3=1; else if type='rectum' then a4=1; else if type='bladder' then a1=1; else a6=1; label type='Type of Cancer' y1='patient 1st hospitialization' x1='control 1st hospitialization' y2='patient survival from untreatability' x2='contorl survival from untreatability' d1='difference 1st hosptialization' d2='difference untreatability' alive='Alive at end of study' ; datalines; stomach 1 f 61 124 264 124 38 0 stomach 2 m 69 42 62 12 18 0 stomach 3 f 62 25 149 19 36 0 stomach 4 f 66 45 18 45 12 0 stomach 5 m 63 412 180 257 64 0 stomach 6 m 79 51 142 23 20 0 stomach 7 m 76 1112 35 128 13 0 stomach 8 m 54 46 299 46 51 0 stomach 9 m 62 103 85 90 10 0 stomach 10 f 69 876 69 876 19 1 stomach 11 m 46 146 361 123 52 0 stomach 12 m 57 340 269 310 28 0 stomach 95 f 59 396 130 359 55 0 bronchus 13 m 74 81 72 74 33 0 bronchus 14 m 74 461 134 423 18 0 bronchus 15 m 66 20 84 16 20 0 bronchus 16 m 52 450 98 450 58 0 bronchus 17 f 48 246 48 87 13 0 bronchus 18 f 64 166 142 115 49 0 bronchus 19 m 70 63 113 50 38 0 bronchus 20 m 77 64 90 50 24 0 bronchus 21 n 71 155 30 113 18 0 bronchus 22 m 70 859 56 857 18 1 bronchus 23 m 39 151 260 38 34 0 bronchus 24 m 70 166 116 156 20 0 bronchus 25 m 70 37 87 27 27 0 bronchus 91 m 55 223 69 218 32 0 bronchus 93 m 74 138 100 138 27 0 bronchus 97 m 69 72 315 39 39 0 bronchus 98 m 73 245 188 231 65 0 colon 28 f 76 248 292 135 18 0 colon 29 f 58 377 492 50 30 0 colon 30 m 49 189 462 189 65 0 colon 31 m 69 1843 235 1267 17 0 colon 32 f 70 180 294 155 57 0 colon 33 f 68 537 144 534 16 0 colon 34 m 50 519 643 502 25 0 colon 35 f 74 455 301 126 21 0 colon 36 m 66 406 148 90 17 0 colon 37 f 76 365 641 365 42 0 colon 38 f 56 942 272 911 40 0 colon 77 m 65 776 198 743 14 1 colon 84 f 74 372 37 366 28 0 colon 90 m 58 163 199 156 31 0 colon 92 f 60 101 154 99 28 0 colon 99 m 77 20 649 20 33 0 colon 100 m 38 283 162 274 80 0 rectum 39 f 56 185 422 62 38 0 rectum 40 f 75 479 82 226 10 0 rectum 41 f 57 875 551 437 62 0 rectum 42 m 56 115 140 85 13 0 rectum 43 m 68 362 106 122 36 0 rectum 44 m 54 244 645 198 80 0 rectum 45 m 59 2175 407 759 64 0 ovary 46 f 49 1234 307 230 64 0 ovary 47 f 68 89 690 88 21 0 ovary 48 f 52 201 285 196 129 0 ovary 49 f 67 356 244 337 58 0 ovary 50 f 56 2970 371 154 66 0 ovary 96 f 51 456 368 91 70 0 breast 51 f 56 1235 796 10 45 0 breast 52 f 57 24 977 24 68 0 breast 53 f 53 1581 1623 580 23 0 breast 54 f 68 1166 555 747 12 0 breast 55 f 68 40 1304 39 52 0 breast 56 f 53 727 1165 87 28 0 breast 57 f 75 3808 675 633 62 0 breast 58 f 68 791 871 251 75 0 breast 59 f 55 1804 916 389 104 0 breast 60 f 43 3460 1311 2270 41 1 breast 61 f 53 719 978 322 64 0 bladder 62 m 93 4288 464 260 29 0 bladder 63 f 70 3658 694 305 22 0 bladder 64 f 77 51 221 37 21 0 bladder 65 f 72 278 490 109 16 0 bladder 66 m 44 548 433 37 11 0 bladder 67 m 64 1607 484 1320 14 1 bladder 68 m 63 1250 152 419 32 1 kidney 71 f 71 205 332 8 91 0 kidney 72 f 63 538 377 96 47 0 kidney 73 f 51 203 147 190 35 0 kidney 74 m 53 296 500 64 34 0 kidney 75 m 57 870 299 260 19 0 kidney 76 m 73 331 585 326 37 0 kidney 78 m 69 1685 1056 46 15 0 kidney 79 m 74 2060 647 2060 44 1 proc sort; by type; proc means n mean; class type; var y1 x1 y2 x2 d1 d2 ; run; proc glm data=vitC; class type ; model d1 d2 = type /solution; * Note: The order of the values in the contrast are alphabetica, so in this case it is bladder bronchus colon kidney rectum stomach; contrast 'bronchus=colon=kidney=rectum=stomach ' type 0 1 0 0 0 -1, type 0 0 1 -1 0 0, type 0 1 -1 -1 0 1, type 0 1 1 1 -4 1; contrast 'bladder vs others' type -5 1 1 1 1 1 ; manova h=type /printh printe; lsmeans type; title 'MANOVA of vitamin C and Cancer'; run; /*************** Simultaneous Confidence Statements ***/ proc glm data=vitC; class type ; model d1 d2 = type /solution; contrast 'bronchus=colon=kidney=rectum=stomach ' type 0 1 0 0 0 -1, type 0 0 1 -1 0 0, type 0 1 -1 -1 0 1, type 0 1 1 1 -4 1; contrast 'bladder vs others' type 5 -1 -1 -1 -1 -1 ; manova h=type m=(1 0) /printh printe; estimate 'b vs o' type 1 -.2 -.2 -.2 -.2 -.2; lsmeans type; title1 'MANOVA of vitamin C and Cancer'; title2 'New Variable: just d1'; run; /****************** MANOVA via IML ***********/ proc iml; /************ Handy modules ***************/ start samplestats(X,Xbar,W,S,R,n); n = nrow(X); one = J(n,1); Xbar = X`*one/n; W = (X - one*Xbar`)` * (X - one*Xbar`); S = W/(n-1); Dsqrt = sqrt(diag(S)); R = inv(Dsqrt)*S*inv(Dsqrt); Finish samplestats; /**********************************************/ /*********** Traditional Method ***************/ use vitC; read all var{d1 d2} where(type='stomach') into Xstomach; read all var{d1 d2} where(type='bronchus') into Xbronchus; read all var{d1 d2} where(type='colon') into Xcolon; read all var{d1 d2} where(type='rectum') into Xrectum; read all var{d1 d2} where(type='bladder') into Xbladder; read all var{d1 d2} where(type='kidney') into Xkidney; run samplestats(Xstomach,Xbars,W,Ss,R,ns); run samplestats(Xbronchus,Xbarb,W,Sb,R,nb); run samplestats(Xcolon,Xbarc,W,Sc,R,nc); run samplestats(Xrectum,Xbarr,W,Sr,R,nr); run samplestats(Xbladder,Xbard,W,Sd,R,nd); run samplestats(Xkidney,Xbark,W,Sk,R,nk); /* Spool */ Wpool = (ns-1)*Ss + (nb-1)*Sb + (nc-1)*Sc + (nr-1)*Sr + (nd-1)*Sd + (nk-1)*Sk; Spool = Wpool/(ns+nb+nc+nr+nd+nk-6); print 'Ellipses using ' Spool; /* MANOVA: Hypothesis Test */ * Have Within groups sscp which i called Wpool ; read all var{d1 d2} into X; run samplestats(X,Xbar,Wtotal,S,R,n); detWpool = det(Wpool); detWtot= det(Wtotal); Wilk = detWpool/detWtot; g = 6; nu_e = n-g; nu_h = g-1; F = ((nu_e-1)/nu_h)*((1 - sqrt(Wilk))/sqrt(Wilk)); df_num = 2*nu_h; df_den = 2*(nu_e-1); pvalue = 1 - cdf('F',F,df_num,df_den); print '1-way MANOVA: ', Wpool, Wtotal, detWpool detWtot , n g nu_e nu_h Wilk, df_num df_den F pvalue; print Xbar; print Spool; /**************************** This is nice to get *********************/ Xbar_groups = Xbard` // Xbarb` // Xbarc` // Xbark` // Xbarr` //Xbars`; Xbar_all = Xbar`; do el = 1 to (g-1); Xbar_all = Xbar_all // Xbar`; end; type = {'Bladder', 'Bronchus', 'Colon', 'Kidney', 'Rectum', 'Stomach' }; Tau = Xbar_groups - Xbar_all; print type Xbar_groups, type Tau; /********* MANOVA as a multivariate General Linear model ********/ print '********* MANOVA as a multivariate General Linear model ********'; print '********** Much simpler, powerful & flexible *******************'; read all var{d1 d2} into X; read all var{ao a1 a2 a3 a4 a5 a6} into A; B = ginv(A`*A)*A`*X; E = X`*X - B`*A`*A*B; run samplestats(X,Xbar,T,S,R,n); detE = det(E); Wilk = detE/det(T); n= nrow(X); p= ncol(X); g= nrow(B)-1; nu_e= n-g; nu_h= g-1; * This is for p=2; F = ((nu_e-1)/nu_h)*((1 - sqrt(Wilk))/sqrt(Wilk)); df_num = 2*nu_h; df_den = 2*(nu_e-1); pvalue = 1 - cdf('F',F,df_num,df_den); type = {'intercept', 'Bladder', 'Bronchus', 'Colon', 'Kidney', 'Rectum', 'Stomach' }; print 'MANOVA Test for no overall effect', type B, n p g nu_e nu_h Wilk, df_num df_den F pvalue, ' '; /* Contrast statement: Bladder vs rest*/ C1={0 -5 1 1 1 1 1}; H1 = (C1*B)`*inv(C1*ginv(A`*A)*C1`) * (C1*B); Wilk1 = detE/det(E+H1); nu_h=nrow(C1); * For nu_h=1; F1=((nu_e+nu_h-p)/p)*((1-Wilk1)/Wilk1); nu_num = p; nu_den = nu_e+nu_h-p; pvalue = 1 -cdf('F',F1,nu_num,nu_den); print 'Testing bladder equals all the others', Wilk1 nu_h nu_e nu_num nu_den F1 pvalue ; /*********** Check that all equal except bladder ********/ C2 = { 0 0 1 -1 0 0 0, 0 0 1 0 -1 0 0, 0 0 1 0 0 -1 0, 0 0 1 0 0 0 -1}; H2 = (C2*B)` * inv(C2*ginv(A`*A)*C2`) * (C2*B); Wilk2 = detE/det(E+H2); nu_h = nrow(C2); nu_num = 2*nu_h; nu_den = 2*(nu_e-1); F2 = ((nu_e-1)/nu_h)*((1-sqrt(Wilk2))/sqrt(Wilk2)); pvalue = 1 - cdf('F',F2,nu_num,nu_den); print 'Testing bronchus = colon = kidney = rectum = stomach', Wilk2 nu_h nu_num nu_den F2 pvalue; /***************** Simultaneous CI-> univariate *****************/ C={0 1 -.2 -.2 -.2 -.2 -.2}; M={1, 0}; diff = C*B*M; se = sqrt( M`*spool*M*C*ginv(A`*A)*C`); F = quantile('F',.95,1,nu_e); lower = diff - sqrt(F)*se; upper = diff + sqrt(F)*se; print diff F se lower upper; /************* Discriminat Functon **************/ proc discrim; class type; var d1 d2; run; proc princomp cov; var d1 d2; run;