'*****************QUANTITATIVE GENETICS PROGRAM********************* 'Written by David Duffy QIMR 1989-92 defdbl a-z color 15,0 MAXTAB=12 ' Maximum number 2x2 tables MAXCAT=15 ' Maximum dimensions for RxC tables DIM num(MAXCAT,MAXCAT),re(MAXCAT),ft(MAXCAT),actual(MAXCAT),expected(MAXCAT) dim xo(2,2,MAXTAB),xn(2,2,MAXTAB),m(MAXTAB),var(MAXTAB) dim flag%(40),stack(40),vars(9),op$(22) dim mon$(12) for i%=1 to 12 read mon$(i%) next i% data "Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec" open "mystat.out" for output as #2 ?#2," +---------------------------------+ " ?#2," | MYSTAT | " ?#2," | Epidemiological statistics | " ?#2," +---------------------------------+ " ?#2," Written by David L Duffy 1989-92 " ?#2," QIMR Australia " ?#2," TurboBasic version " ?#2, ?#2," MYSTAT run at ";time$;" on ";mid$(date$,4,3);mon$(val(left$(date$,2)));_ right$(date$,5) ?#2, 100 cls:?#2,: ?" Welcome to my statistics program" 110 ?:?" We currently have the following on the menu:" 120 ?" (1) Chi-square goodness-of-fit ?" (2) Contingency Chi-square (incl 2xK analysis) ?" (3) Chi-square K proportions (2xK analysis) ?" (4) Symmetry model for square contingency table ?" (5) Evaluate central Chi-square ?" (6) Byar test for SMRs ?" (7) z-test for incidences ?" (8) z-test for proportions ?" (9) odds ratio ?" (10) Fisher-Irwin Exact test & exact CI for OR ?" (11) Matched analysis OR and exact CI ?" (12) Tetrachoric correlation ?" (13) Cohen's Kappa ?" (14) Exact CI for proportion ?" (15) Multiple 2x2 tables ?" (16) Calculator ? 200 INPUT"Enter the number of your choice: ", choice 210 ON choice GOSUB 11000,11500,11700,18000,10900,1200,1000,_ 2000,1800,13000,22000,14000,16000,17000,20000,21000 220 ? 250 ?"Another one? [y/N/a (same again)]" 260 a$=INKEY$:IF a$="" THEN 260 265 if a$="a" or a$="A" then 210 if a$>"0" and a$<":" then choice=val(a$):goto 210 270 IF a$="Y" OR a$="y" THEN 120 ELSE END 1000 ' z-test for person-time data from Rothman cls:?#2, ?"----------------------------------------------------" ?#2,"----------------------------------------------------" ?" Z-test for two incidence rates" ?#2," Z-test for two incidence rates" ?"----------------------------------------------------" ?#2,"----------------------------------------------------" ?"Enter number of positives in group one : ";:INPUT, pos1 ?#2,"Enter number of positives in group one : ";pos1 ?"Enter person-time observed for group one : ";:INPUT, n1 ?#2,"Enter person-time observed for group one : ";n1 ?"Enter number of positives in group two : ";:INPUT, pos2 ?#2,"Enter number of positives in group two : ";pos2 ?"Enter person-time observed for group two : ";:INPUT, n2 ?#2,"Enter person-time observed for group two : ";n2 sed=sqr((pos1+pos2)*n1*n2)/(n1+n2) z=(pos1-n1*(pos1+pos2)/(n1+n2))/sed ? ?#2, ?"Incidence rate in group one= ";:? USING "#.#####";pos1/n1 ?#2,"Incidence rate in group one= ";:?#2, USING "#.#####";pos1/n1 ?"Incidence rate in group two= ";:? USING "#.#####";pos2/n2 ?#2,"Incidence rate in group two= ";:?#2, USING "#.#####";pos2/n2 ? ?#2, rd=pos1/n1-pos2/n2 sed=sqr(pos1/n1/n1+pos2/n2/n2) w=pos1*n2/pos2/n1 logse=sqr(1/pos1+1/pos2) ?using"Rate Difference= +#.#####";rd; ?#2,using"Rate Difference= +#.#####";rd; ?"; 95% CI= ";:?USING " +#.###### to +#.######";rd-1.96*sed;rd+1.96*sed ?#2,"; 95% CI= ";:?#2,USING " +#.###### to +#.######";rd-1.96*sed;rd+1.96*sed ?using" Rate Ratio= ###.##";w; ?#2,using" Rate Ratio= ###.##";w; ?"; 95% CI= ";:? USING "###.##";exp(log(w)+1.96*logse); ?#2,"; 95% CI= ";:?#2, USING "###.##";exp(log(w)+1.96*logse); ?" to ";:? USING"###.##";exp(log(w)-1.96*logse) ?#2," to ";:?#2, USING"###.##";exp(log(w)-1.96*logse) ? ?#2, ?using" z= ##.##"; abs(z); ?#2,using" z= ##.##"; abs(z); call zprob(z) RETURN ' 1200 'Byars test for approximating poisson distribution for SMRs ?"----------------------------------------------------" ?#2,"----------------------------------------------------" ?" Byars Test for an SMR " ?#2," Byars Test for an SMR " ?"----------------------------------------------------" ?#2,"----------------------------------------------------" ? ?#2, ?"Enter OBSERVED number of events: ";:INPUT, actual(1) ?#2,"Enter OBSERVED number of events: ";actual(1) ?"Enter EXPECTED number of events: ";:INPUT, expected(1) ?#2,"Enter EXPECTED number of events: ";expected(1) smr=actual(1)/expected(1) if actual(1)=0 then w=p1*(1-p2)/p2/(1-p1) else w=1e12 logse=1/(pos1+0.5)+1/(n1-pos1+0.5)+1/(pos2+0.5)+1/(n2-pos2+0.5) logse=SQR(logse) relrisk=pos1*n2/pos2/n1 relrse=sqr((n1-pos1)/n1/pos1+(n1-pos2)/n2/pos2) ?:?#2, ?using" Proportion positive in group one= #.##";p1 ?#2,using" Proportion positive in group one= #.##";p1 ?using" Proportion positive in group two= #.##";p2 ?#2,using" Proportion positive in group two= #.##";p2 ? ?#2, ?" Type of Measure Estimate Standard 95% CL Cornfield 95% CL " ?#2," Type of Measure Estimate Standard 95% CL Cornfield 95% CL " ?" --------------- -------- ---------------- ------------------ " ?#2," --------------- -------- ---------------- ------------------ " ?using " Rate Difference +#.##";p1-p2; ?#2,using " Rate Difference +#.##";p1-p2; ?using " +#.## to -#.##";p1-p2-1.96*sed;p1-p2+1.96*sed; ?#2,using " +#.## to -#.##";p1-p2-1.96*sed;p1-p2+1.96*sed; call corny(pos1,n1-pos1,pos2,n2-pos2,ul,ll,3) ?using " +#.## to +#.##";ll;ul ?#2,using " +#.## to +#.##";ll;ul ?using " Odds Ratio ###.##";w; ?#2,using " Odds Ratio ###.##";w; if w<>0 then ll=EXP(LOG(w)-1.96*logse) else ll=0 if w<>0 then ul=EXP(LOG(w)+1.96*logse) else ul=exp(1e-12+1.96*logse) ?using " ###.## to ###.##";ll;ul; ?#2,using " ###.## to ###.##";ll;ul; call corny(pos1,n1-pos1,pos2,n2-pos2,ul,ll,1) ?using " ###.## to ###.##";ll;ul ?#2,using " ###.## to ###.##";ll;ul ll=exp(log(relrisk)-1.96*relrse) ul=exp(log(relrisk)+1.96*relrse) ?using " Relative Risk ###.##";relrisk; ?#2,using " Relative Risk ###.##";relrisk; ?using " ###.## to ###.##";ll;ul; ?#2,using " ###.## to ###.##";ll;ul; call corny(pos1,n1-pos1,pos2,n2-pos2,ul,ll,2) ?using " ###.## to ###.##";ll;ul ?#2,using " ###.## to ###.##";ll;ul ?using " Phi Coefficient +#.##";phiz; ?#2,using " Phi Coefficent +#.##";phiz; call corny(pos1,n1-pos1,pos2,n2-pos2,ul,ll,4) ?using " +#.## to +#.##";ll;ul ?#2,using " +#.## to +#.##";ll;ul ? ?#2, ? ?#2, ?using" Binomial z test= ##.##"; abs(z); ?#2,using" Binomial z test= ##.##"; abs(z); call zprob(z) 2160 RETURN 10900 '********************evaluate chi-square******************* ?#2, input "Enter Chi-square value : ",chisq ?#2, "Chi-square value : ";chisq input "Enter degrees of freedom: ",df ?#2, "degrees of freedom: ";df ?:?using"AIC=+####.#";chisq-2*df; ?#2,:?#2,using"AIC=+####.#";chisq-2*df; call chiprob(chisq,df) return 11000 '********************goodness-of-fit*********************** ?"----------------------------------------------" ?#2,"----------------------------------------------" ?" Goodness-of-fit Chi-square test " ?#2," Goodness-of-fit Chi-square test " ?"----------------------------------------------" ?#2,"----------------------------------------------" 11010 INPUT "How many categories to be tested: ",cats ?#2, "How many categories to be tested: ";cats 11025 total=0 11030 FOR c=1 TO cats 11040 ?"Enter expected value no. ";c;:INPUT ": ",expected(c) ?#2,"Enter expected value no. ";c;": ";expected(c) 11050 ?"Enter actual value no. ";c;:INPUT ": ",actual(c) ?#2,"Enter actual value no. ";c;": ";actual(c) 11060 total=total+actual(c) 11070 NEXT c 11110 chisq=0 11115 cls:?#2,: ?"------------- Chi-squared goodness of fit ------------------" cls:?#2,: ?#2,"------------- Chi-squared goodness of fit ------------------" 11117 ?" OBS EXP DEVIATE (F-T) " ?#2," OBS EXP DEVIATE (F-T) " 11120 FOR c=1 TO cats 11130 a=actual(c)-expected(c) 11135 ft=sqr(actual(c))+sqr(actual(c)+1)-sqr(4*expected(c)+1) 11140 chisq=chisq+a*a/expected(c) 11150 ? USING "####.## "; actual(c); expected(c) ; ft ?#2, USING "####.## "; actual(c); expected(c) ; ft 11170 NEXT c 11180 ?"------------------------------------------------------------" ?#2,"------------------------------------------------------------" 11185 df=cats-1 11190 ?"G.O.F. X";chr$(253);"= "; ?#2,"G.O.F. X";chr$(253);"= "; 11195 ? USING "###.##";chisq;:?" (df= ";df;")"; ?#2, USING "###.##";chisq;:?#2," (df= ";df;")"; 11210 call chiprob(chisq,df) 11215 RETURN 11500 '**********************contingency chi sq******************** 11510 cls:?#2, ?"-------------------------------------------------" ?#2,"-------------------------------------------------" ?" CONTINGENCY CHI-SQUARE TEST " ?#2," CONTINGENCY CHI-SQUARE TEST " ?"-------------------------------------------------" ?#2,"-------------------------------------------------" INPUT "How many categories for the row variable : ",numa ?#2, " No. of categories for the row variable : ";numa INPUT "How many categories for the column variable: ",numb ?#2, " No. of categories for the column variable: ";numb cls:?#2, for r=1 to numa:for c=1 to numb ?" ";:color 0,7:? r;c;:color 15,0:?" "; ?#2," ";r;c;": "; input;, num(r,c) ?#2,num(r,c) next c ?:?#2,:next r if numa=2 then for r=1 to numa: for c=r to numb swap num(r,c), num(c,r) next:next swap numa, numb end if goto 12000 11700 '********************** enter K proportions ***************** ?"-------------------------------------------------" ?#2,"-------------------------------------------------" ?" CONTINGENCY CHI-SQUARE TEST K PROPORTIONS " ?#2," CONTINGENCY CHI-SQUARE TEST K PROPORTIONS " ?"-------------------------------------------------" ?#2,"-------------------------------------------------" input "How many proportions: ",numa numb=2 ? ?#2, for r=1 to numa ?"Group ";:?using "## ";r; ?#2,"Group ";:?#2,using "## ";r; ?" ";:color 0,7:? " events";:color 15,0:?" ";:input;, num(r,1) ?#2," events: ";num(r,1); ?" ";:color 0,7:? " denom ";:color 15,0:?" ";:input, num(r,2) ?#2," denom : ";num(r,2) num(r,2)=num(r,2)-num(r,1) next r 12000 cls:?#2, ?"-------------------- CONTINGENCY CHI-SQ ---------------------" ?#2,"-------------------- CONTINGENCY CHI-SQ ---------------------" if numb=2 then ?"Gp Events Non-Events Rate Odds Ratio Cornfield 95% CL" ?#2,"Gp Events Non-Events Rate Odds Ratio Cornfield 95% CL" end if FOR r=1 TO numa:suma(r)=0:NEXT r FOR c=1 TO numb:sumb(c)=0:NEXT c sumtot=0:sr=0:sc=0 FOR r=1 TO numa:FOR c=1 TO numb mr=mr+(r-1)*num(r,c) mc=mc+(c-1)*num(r,c) suma(r)=suma(r)+num(r,c) sumb(c)=sumb(c)+num(r,c) sumtot=sumtot+num(r,c) NEXT:NEXT mr=mr/sumtot mc=mc/sumtot chisq=0:gsq=0:ftsq=0:cov=0 FOR r=1 TO numa:FOR c=1 TO numb expfreq=suma(r)*sumb(c)/sumtot cov=cov+num(r,c)*(r-1-mr)*(c-1-mc) sr=sr+num(r,c)*(r-1-mr)*(r-1-mr) sc=sc+num(r,c)*(c-1-mc)*(c-1-mc) ft(c)=sqr(num(r,c))+sqr(num(r,c)+1)-sqr(4*expfreq+1) chisq=chisq+(num(r,c)-expfreq)^2/expfreq if num(r,c)<>0 then gsq=gsq+2*num(r,c)*log(num(r,c)/expfreq) ftsq=ftsq+ft(c)*ft(c) re(c)=expfreq next c if numb=2 then ? "#";:?using "##";r; ?#2, "#";:?#2,using "##";r; else ? "Obs "; ?#2, "Obs "; end if for c=1 to numb ? using " ###### ";num(r,c); ?#2, using " ###### ";num(r,c); next c if numb=2 and r=1 then ? using " #.##### ###.# ";num(r,1)/(num(r,1)+num(r,2)),1 ?#2, using " #.##### ###.# ";num(r,1)/(num(r,1)+num(r,2)),1 elseif numb=2 and r>1 then if num(r,2)=0 or num(1,1)=0 then oddsr=999.9 else oddsr=num(r,1)*num(1,2)/num(r,2)/num(1,1) end if call corny(num(r,1),num(r,2),num(1,1),num(1,2),ul,ll,1) ? using " #.##### ###.# ###.# to ###.# "; _ num(r,1)/(num(r,1)+num(r,2)),oddsr,ll,ul ?#2, using " #.##### ###.# ###.# to ###.# "; _ num(r,1)/(num(r,1)+num(r,2)),oddsr,ll,ul else ? ?#2, end if if numb>2 then ? "Exp "; ?#2, "Exp "; for c=1 to numb ? using "{####.##}";re(c); ?#2, using "{####.##}";re(c); next c ? ?#2, ? "Dev "; ?#2, "Dev "; for c=1 to numb ? using "[ +##.##]";ft(c); ?#2, using "[ +##.##]";ft(c); next c ? ?#2, end if next r mh=(sumtot-1)*cov*cov/sr/sc df=(numa-1)*(numb-1) ?"-------------------------------------------------------------" ?#2,"-------------------------------------------------------------" ? " X";chr$(253);" = ";:?using "###.##";chisq; ?#2, " X";chr$(253);" = ";:?#2,using "###.##";chisq; ? " (";df;"df)"; ?#2, " (";df;"df)"; call chiprob(chisq,df) ? " G";chr$(253);" = ";:?using "###.##";gsq;:chisq=gsq ?#2, " G";chr$(253);" = ";:?#2,using "###.##";gsq;:chisq=gsq call chiprob(chisq,df) ? " Z";chr$(253);" = ";:?using "###.##";ftsq;:chisq=ftsq ?#2, " Z";chr$(253);" = ";:?#2,using "###.##";ftsq;:chisq=ftsq call chiprob(chisq,df) ? "Trend X";chr$(253);" = ";:?using "###.##";mh;:chisq=mh ?#2, "Trend X";chr$(253);" = ";:?#2,using "###.##";mh;:chisq=mh ? " ( 1 df)"; ?#2, " ( 1 df)"; call chiprob(chisq,1) return 13000 '*******************Fisher-Irwin Exact Test******************* 'Exact two sided Confidence Limits for Odds Ratio (2x2 table) ' 'Algorithm AS115 Applied Statistics 26:214-220 by Baptista J, Pike MC. ' 'ported to TurboBasic from original FORTRAN by DL Duffy ' defdbl a-h,k-z defint i,j ' 'Entry routine ' cls:?#2, ?"-------------------------" ?#2,"-------------------------" ?" Fisher Exact Test 2x2" ?#2," Fisher Exact Test 2x2" ?"-------------------------" ?#2,"-------------------------" call twobytwo(a,b,c,d) itab(1)=a itab(2)=b itab(3)=c itab(4)=d ' 'Initialize ' n12max=0 'length of vectors of log factorials maxlf=500 'max limit for a+b or c+d big=99999 alpha=0.05 ladz=0 lbcz=0 ' dim lf(maxlf) 'workspace log factorials dim p(maxlf) 'workspace vector of probs for x given OR dim itab(4),jtab(4) ' n1=itab(1)+itab(2) n2=itab(3)+itab(4) m=itab(1)+itab(3) jtab(1)=itab(2) jtab(2)=itab(1) jtab(3)=itab(4) jtab(4)=itab(3) ' 'calculates approx limits ' v=exp(sqr(1/(a+0.5)+1/(b+0.5)+1/(c+0.5)+1/(d+0.5))*1.96) ' if itab(1)=0 or itab(4)=0 then ladz=1 if itab(2)=0 or itab(3)=0 then lbcz=0 if ladz=1 then oddsr=0 all=0 aul=1 elseif lbcz=1 then oddsr=big aul=big all=1 else oddsr=a*d/b/c all=oddsr/v aul=oddsr*v end if ' if n1>n2 then nmax=n1+1 else nmax=n2+1 if nmax>maxlf then ?" TOTALS TOO LARGE, use asymptotic test ":erase lf,p,itab,jtab:return if nmax>maxlf then ?#2," TOTALS TOO LARGE, use asymptotic test ":erase lf,p,itab,jtab:return if nmax=2 then i2=n12max else i2=2 e=i2-2 ' for i=i2 to nmax e=e+1 lf(i)=log(e)+lf(i-1) next i ' n12max=nmax ' 13060 start=1 if ladz=1 then 13110 if lbcz=1 then eul=big if ladz=1 then 13120 ' 'lower limit ' call itera(itab(),alpha,maxlf,lf(),p(),all,ell) if lbcz=1 then 13120 ' 13090 start=1/aul ' 'upper limit ' call itera(jtab(),alpha,maxlf,lf(),p(),start,eul) eul=1/eul goto 13120 ' 13110 ell=0 goto 13090 ' 'compute exact upper and lower tail areas for OR=1 ' 13120 call proba(1,n1,n2,m,maxlf,lf(),p(),minj,maxj) j=itab(1)+1 if oddsr>1 then 13140 elp=0 ' for i=minj to j elp=elp+p(i) next i ' eup=1-elp+p(j) goto 13200 ' 13140 eup=0 ' for i=j to maxj eup=eup+p(i) next i ' elp=1-eup+p(j) 13200' ' 'Output results ' ? ?#2, ?using" Odds ratio= ###.####";oddsr ?#2,using" Odds ratio= ###.####";oddsr ?using"Exact 95% lower limit= ###.####";ell ?#2,using"Exact 95% lower limit= ###.####";ell ?using"Exact 95% upper limit= ###.####";eul ?#2,using"Exact 95% upper limit= ###.####";eul ? ?#2, ?"EXACT TEST (ie for OR=1)" ?#2,"EXACT TEST (ie for OR=1)" ?using"Lower tail area = ###.####";elp ?#2,using"Lower tail area = ###.####";elp ?using"Upper tail area = ###.####";eup ?#2,using"Upper tail area = ###.####";eup ' erase lf,p,itab,jtab return ' 'Subroutines ' sub inout(in,ia,maxlf,p(1),minj,maxj,alpha) ' 'gives 1 if ia in given confidence interval ' local j,j1,j2,sum,pmin in=1 j1=minj j2=maxj j=ia+1 sum=0 ' ' adds terms from either tail in ascending order of size ' 13520 if p(j1)=alpha then 13550 if pmin=p(j2) then 13530 j1=j1+1 if j1>j then 13540 goto 13520 ' 13530 j2=j2-1 if j2>=j then 13520 13540 in=0 13550 end sub ' ' sub itera(itab(1),alpha,maxlf,lf(1),p(1),start,rho) ' 'returns rho as lower limit, start=start value ' local psi,eps eps=0.000001 n1=itab(1)+itab(2) n2=itab(3)+itab(4) m=itab(1)+itab(3) psi=start rho=0 kl=0 ' 13310 call proba(psi,n1,n2,m,maxlf,lf(),p(),minj,maxj) ' call inout(in,itab(1),maxlf,p(),minj,maxj,alpha) ' 'kl=1 when rho and opsi span correct value ' if kl=1 then 13340 if in=1 then 13320 rho=psi psi=psi*1.1 goto 13310 13320 kl=1 opsi=psi 13330 psi=(rho+opsi)*0.5 ' 'new estimate is midpoint of spanning interval ' goto 13310 13340 if in=1 then opsi=psi else rho=psi if abs(rho/opsi-1)0 then minj=m-n2+1 else minj=1 if m>n1 then maxj=n1+1 else maxj=m+1 alt=log(psi) if abs(psi-1)maxj then maxpj=maxj goto 13220 ' 13210 maxpj=fix(m/(n1+n2)*n1+1.5) 13220 j=maxpj n1j2=n1-j+2 mj2=m-j+2 n2mj=n2-m+j ' 'compute pmax, log of the unscaled unconditional 'prob at maxpj ' pmax=lf(n1+1)-lf(j)-lf(n1j2)+lf(n2+1)-lf(mj2)-lf(n2mj)+(j-1)*alt p(j)=1 ps=1 temp=lf(n1+1)+lf(n2+1)-pmax maxpj=maxpj-1 if maxpjeps then 13240 i=l+1 if i>maxpj then 13250 for ll=i to maxpj j=j-1 p(j)=0 next ll goto 13250 13240 ' ' next l ' 13250 maxpj=maxpj+2 if maxpj>maxj then 13280 n1j2=n1-maxpj+2 mj2=m-maxpj+2 n2mj=n2-m+maxpj for j=maxpj to maxj p(j)=temp-lf(j)-lf(n1j2)-lf(mj2)-lf(n2mj)+(j-1)*alt p(j)=exp(p(j)) ps=ps+p(j) n1j2=n1j2-1 mj2=mj2-1 n2mj=n2mj+1 if p(j)>eps then 13270 i=j+1 if i>maxj then 13280 for l=i to maxj p(l)=0 next l goto 13280 13270 ' next j ' 'compute conditional probabilities by scaling ' 13280 for j=minj to maxj p(j)=p(j)/ps next j end sub ' 'Kirk's algorithm for tetrachoric correlation ' ported from original FORTRAN ' Psychometrika 1979;44:461 ' 14000 '*********************TETRACHORIC CORRELATION********************* ?"-------------------------------------------------" ?#2,"-------------------------------------------------" ?" TETRACHORIC CORRELATION " ?#2," TETRACHORIC CORRELATION " ?"-------------------------------------------------" ?#2,"-------------------------------------------------" call twobytwo(a,b,c,d) a1=a 14102 if d0.8) THEN est=SGN(a)*0.8 14670 'converg c for r 14680 eps=0.0001 14690 GOTO 14760 14700 IF x>0.5 THEN 15020 14710 'Hastings approximation 14720 e=SQR(-2*LOG(x)) 14730 e1=(0.010328*e+0.802853)*e+2.515517 14740 e2=((0.001308*e+0.189269)*e+1.432788)*e+1 14750 est=e-e1/e2 14760 old=est 14770 'main loop 14780 FOR i=1 TO 20 14789 x=old 14790 xx=old*old 14800 IF j=3 GOTO 14860 14810 'h,k integrand evaluation for gaussian quadrature 14820 f=w1*(FNone(c1sq)+FNone(c8sq))+w2*(FNone(c2sq)+FNone(c7sq))+w3*(FNone(c3sq)+FNone(c6sq)) 14830 f=f+w4*(FNone(c4sq)+FNone(c5sq)):funcnum=old*f 14840 funcnew=old-(funcnum-con)/FNone(1) 14850 GOTO 14930 14860 IF ABS(old)>1 THEN 14980 14870 'r integrand evaluator 14880 g=w1*(FNtwo(c1,c1sq)+FNtwo(c8,c8sq))+w2*(FNtwo(c2,c2sq)+FNtwo(c7,c7sq)) 14890 g=g+w3*(FNtwo(c3,c3sq)+FNtwo(c6,c6sq))+w4*(FNtwo(c4,c4sq)+FNtwo(c5,c5sq)) 14900 funcnum=old*g 14910 funcnew=old-(funcnum-con)/FNtwo(1,1) 14920 'test for convergence 14930 IF ABS(old-funcnew)>eps THEN 14950 14940 ON j GOTO 14470,14540,15040 14950 old=funcnew 14955 NEXT i 14960 r=3 14970 GOTO 15020 14980 IF r=2 THEN 15020 14990 r=2 15000 old=SGN(a)*0.97 15010 GOTO 14780 15020 r=4 15030 GOTO 15045 15040 r=funcnew 15045 IF r>1 THEN ?"Marginal proportion >0.5;try other diagonal":RETURN IF r>1 THEN ?#2,"Marginal proportion >0.5;try other diagonal":RETURN 15050 ? using"tetrachoric r= +#.###";r ?#2, using"tetrachoric r= +#.###";r 15052 ' Asymptotic SE 15053 phi=exp(-(h*h-2*r*h*fk+fk*fk)/2/(1-r*r))/6.2831853/sqr(1-r*r) 15054 rec=1/sqr(1/a1+1/b+1/c+1/d) 15055 se=rec/n/phi 15056 ?"asymptotic standard error (r)="; ?#2,"asymptotic standard error (r)="; 15057 ?using "#.###";se ?#2,using "#.###";se 15060 RETURN 16000 ' Fleiss, Cohen & Everitt asymptotic estimator for variance ' of kappa supplemented by bootstrap. cls:?#2, ?" COHEN KAPPA FOR TWO VARIABLES " ?#2," COHEN KAPPA FOR TWO VARIABLES " INPUT "How many categories for variables: ",n ?#2, " No. of categories for variables: ";n cls:?#2, tot=0:var=0:mk=0 for r=1 to n:for c=1 to n ?" ";:color 0,7:? r;c;:color 15,0:?" "; input;, num(r,c) ?#2," ";r;c;" ";num(r,c) tot=tot+num(r,c) next c ?:?#2,:next r call korig(num(),tot,n,kap,ase) kappa=kap if ase>0 then ase=sqr(ase) else ase=0 for r=1 to n for c=1 to n num(r,c)=num(r,c)-1 call kappa(num(),n,kap) k(r,c)=tot*kappa-(tot-1)*kap num(r,c)=num(r,c)+1 mk=mk+num(r,c)*k(r,c) next:next mk=mk/tot for r=1 to n for c=1 to n var=var+num(r,c)*(k(r,c)-mk)*(k(r,c)-mk) next:next bse=sqr(var/tot/(tot-1)) ? ?#2, ?using" Cohen Kappa = +#.##";kappa ?#2,using" Cohen Kappa = +#.##";kappa ?using" Jackknife K = +#.##";mk ?#2,using" Jackknife K = +#.##";mk ?using" Asymp SE = #.##";ase ?#2,using" Asymp SE = #.##";ase ?using" Jack' SE = #.##";bse ?#2,using" Jack' SE = #.##";bse ?using" Asym 95% CI = +#.## to +#.##";kappa-1.96*ase;kappa+1.96*ase ?#2,using" Asym 95% CI = +#.## to +#.##";mkappa-1.96*ase;kappa+1.96*ase ?using" Jack 95% CI = +#.## to +#.##";mk-1.96*bse;mk+1.96*bse ?#2,using" Jack 95% CI = +#.## to +#.##";mk-1.96*bse;mk+1.96*bse return ' sub kappa(num(2),n,kap) local r,c,t1,t2,tot,ro(),co() shared MAXCAT dim ro(MAXCAT),co(MAXCAT) t1=0:t2=0 for r=1 to n for c=1 to n ro(c)=ro(c)+num(c,r) co(c)=co(c)+num(r,c) tot=tot+num(r,c) next:next for r=1 TO n t1=t1+num(r,r)/tot t2=t2+co(r)*ro(r)/tot/tot next r kap=(t1-t2)/(1-t2) end sub ' sub korig(num(2),tot,n,kap,var) local p(),ro(),co() shared MAXCAT dim p(MAXCAT,MAXCAT),ro(MAXCAT),co(MAXCAT) t1=0:t2=0:t3=0:t4=0 for r=1 to n for c=1 to n p(r,c)=num(r,c)/tot next:next for r=1 to n for c=1 to n ro(c)=ro(c)+p(c,r) co(c)=co(c)+p(r,c) next:next for r=1 TO n t1=t1+p(r,r) t2=t2+co(r)*ro(r) t3=t3+co(r)*ro(r)*(co(r)+ro(r)) next r kap=(t1-t2)/(1-t2) var=(t2+t2*t2-t3)/tot/(1-t2)/(1-t2) end sub ' ' 17000 ' Pearson-Clopper 95% CL for proportions ?"------------------------------------------" ?#2,"------------------------------------------" ?" Pearson-Clopper 95% CL for proportions" ?#2," Pearson-Clopper 95% CL for proportions" ?"------------------------------------------" ?#2,"------------------------------------------" input "numerator = ",num input "denominator= ",den ?#2, "numerator = ";num ?#2, "denominator= ";den prop=num/den call ci95(num,den,upp95,low95) ? ?#2, ?"Proportion = "; ?#2,"Proportion = "; ?using "#.#######";prop ?#2,using "#.#######";prop ?"Exact 95%CI= "; ?#2,"Exact 95%CI= "; ?using "#.#######"; low95; ?#2,using "#.#######"; low95; ?" to "; ?#2," to "; ?using "#.#######"; upp95 ?#2,using "#.#######"; upp95 return ' ' incomplete beta function ' def FNbetai(a,b,x) bt=exp(FNdlogam(a+b)-FNdlogam(a)-FNdlogam(b)+a*log(x)+b*log(1-x)) if (x<(a+1.0)/(a+b+2.0)) then FNbetai=bt*FNbetacf(a,b,x)/a else FNbetai=1.0-bt*FNbetacf(b,a,1.0-x)/b end if end def ' ' ACM 291 Comm ACM 1966, 9:684 ' def FNdlogam(x) local a1#,a2#,a3#,a4#,a5#,half#,zero#,one#,seven#,y#,f#,z# a1#=0.918938533204673 : a2#=0.000595238095238 a3#=0.000793650793651 : a4#=0.002777777777778 : a5#=0.083333333333333 half#=0.5d0 : zero#=0.0d0 : one#=1.0d0 : seven#=7.0d0 ' FNdlogam=zero# y#=x f#=zero# if (y#ITMAX if m#>ITMAX then print 'ERROR in BETACF - A or B too big, or maxit too small' FNbetacf=az# end def ' ' ' 95% confidence interval for proportion calling betacf ' sub ci95(num,den,upper,lower) ' lim=0.025d0 prop=num/den if num=0 then lower=0.00 upper=1-exp(log(lim)/den) elseif num=den then upper=1.00 lower=exp(log(lim)/den) else ' ' main option ' ste=prop/2 dir=-1 lim=0.025 call it(den-num,den,lim,ste,dir,upper) upper=1.0d0-upper ' ste=prop/2 dir=-1 lim=0.025 call it(num,den,lim,ste,dir,lower) end if end sub ' ' iterate for bound ' sub it(num,den,lim,ste,dir,trial1) ' iter=0 ptail=0 trial1=num/den while (abs(lim-ptail)>0.0000001 and iter<500) iter=iter+1 trial2=trial1+dir*ste if (trial2<0.0000001d0) then trial2=0.0000001d0 if (trial2>0.9999999d0) then trial2=0.9999999d0 ptail=FNbetai(num,den-num+1,trial2) if sgn(dir*(lim-ptail))=1 then trial1=trial2 else ste=ste/2 end if wend end sub 18000 cls:?#2, ?" TESTING OF SYMMETRY MODEL FOR SQUARE CONTINGENCY TABLE" ?#2," TESTING OF SYMMETRY MODEL FOR SQUARE CONTINGENCY TABLE" INPUT "How many categories for variables: ",n ?#2," No. of categories for variables: ";n cls:?#2, for r=1 to n:for c=1 to n ?" ";:color 0,7:? r;c;:color 15,0:?" "; input;, num(r,c) next c ?:next r cls:?#2, ?"---------------------- SYMMETRY MODEL -------------------------" ?#2,"---------------------- SYMMETRY MODEL -------------------------" chisq=0:gsq=0:ftsq=0:upptri=0:lowtri=0 FOR r=1 TO n:FOR c=1 TO n if r<>c then expfreq=(num(r,c)+num(c,r))/2 else expfreq=num(r,c) end if if c>r then upptri=upptri+num(r,c) if r>c then lowtri=lowtri+num(r,c) ft(c)=sqr(num(r,c))+sqr(num(r,c)+1)-sqr(4*expfreq+1) if expfreq>0 then chisq=chisq+(num(r,c)-expfreq)^2/expfreq if num(r,c)>0 then gsq=gsq+2*num(r,c)*log(num(r,c)/expfreq) ftsq=ftsq+ft(c)*ft(c) re(c)=expfreq next c ? "Obs ";:for c=1 to n:? using " #### ";num(r,c);:next c:? ?#2, "Obs ";:for c=1 to n:?#2, using " #### ";num(r,c);:next c:?#2, ? "Exp ";:for c=1 to n:? using "{####.##}";re(c);:next c:? ?#2, "Exp ";:for c=1 to n:?#2, using "{####.##}";re(c);:next c:?#2, ? "Dev ";:for c=1 to n:? using "[ +##.##]";ft(c);:next c:? ?#2, "Dev ";:for c=1 to n:?#2, using "[ +##.##]";ft(c);:next c:?#2, next r df=n*(n-1)/2 ?"---------------------------------------------------------------" ?#2,"---------------------------------------------------------------" ? "Symmetry X";chr$(253);"= ";:?using "###.##";chisq; ?#2, "Symmetry X";chr$(253);"= ";:?#2,using "###.##";chisq; ? " (";df;"df)"; ?#2, " (";df;"df)"; call chiprob(chisq,df) ? "Symmetry G";chr$(253);"= ";:?using "###.##";gsq;:chisq=gsq ?#2, "Symmetry G";chr$(253);"= ";:?#2,using "###.##";gsq;:chisq=gsq call chiprob(chisq,df) ? "Symmetry Z";chr$(253);"= ";:?using "###.##";ftsq;:chisq=ftsq ?#2, "Symmetry Z";chr$(253);"= ";:?#2,using "###.##";ftsq;:chisq=ftsq call chiprob(chisq,df) z=(upptri-lowtri)/sqr(upptri+lowtri) ?using " McNemar Z = +##.##";z; ?#2,using " McNemar Z = +##.##";z; call zprob(z) call Quasisym(n,num()) return ' ' Mantel-Haenzel and ML procedures for 2x2xK tables ' 20000 cls:?#2, ?"------------------------------------------" ?#2,"------------------------------------------" ?" Combination multiple 2x2 tables ?#2," Combination multiple 2x2 tables ?"------------------------------------------" ?#2,"------------------------------------------" tot=0 : chisq=0 mhoddsr=0: mhnum=0: mhden=0: chimh=0: var=0 sumpr=0: sumpsqr=0: sumqs=0: sumr=0: sums=0 ? ?" D+ D- ?" ------------------- ?" E+ a b ?" E- c d ? input "How many tables: ",ns ? ?#2, " No. of tables : ";ns for i=1 to ns ?using"Table ## ";i; ?#2,using"Table ## ";i; ?" ";:color 0,7:? " a ";:color 15,0:?" ";:input;, xo(1,1,i) ?#2," a: ";xo(1,1,i); ?" ";:color 0,7:? " b ";:color 15,0:?" ";:input;, xo(1,2,i) ?#2," b: ";xo(1,2,i); ?" ";:color 0,7:? " c ";:color 15,0:?" ";:input;, xo(2,1,i) ?#2," c: ";xo(2,1,i); ?" ";:color 0,7:? " d ";:color 15,0:?" ";:input, xo(2,2,i) ?#2," d: ";xo(2,2,i) for j=1 to 2: for k=1 to 2 xn(j,k,i)=xo(j,k,i)+0.5 tot=tot+xo(j,k,i) next:next:next for i=1 to ns t=xo(1,1,i)+xo(2,2,i)+xo(1,2,i)+xo(2,1,i) m(i)=xo(1,1,i)-(xo(1,1,i)+xo(2,1,i))*(xo(1,1,i)+xo(1,2,i))/t var(i)=(xo(1,1,i)+xo(2,1,i))*(xo(1,1,i)+xo(1,2,i))*_ (xo(2,1,i)+xo(2,2,i))*(xo(1,2,i)+xo(2,2,i))/t/t/(t-1) sumpr=sumpr+(xo(1,1,i)+xo(2,2,i))*xo(1,1,i)*xo(2,2,i)/t/t sumpsqr=sumpsqr+(xo(1,1,i)*xo(1,2,i)*xo(2,1,i)+xo(1,2,i)*xo(2,1,i)*xo(2,2,i)_ +xo(1,1,i)*xo(1,2,i)*xo(2,2,i)+xo(1,1,i)*xo(2,1,i)*xo(2,2,i))/t/t sumqs=sumqs+(xo(1,2,i)+xo(2,1,i))*xo(1,2,i)*xo(2,1,i)/t/t sumr=sumr+xo(1,1,i)*xo(2,2,i)/t sums=sums+xo(1,2,i)*xo(2,1,i)/t chimh=chimh+m(i) var=var+var(i) next i mhoddsr=sumr/sums call itprop(mloddsr,ns,tot,xo(),xn(),chisq) call jack(mloddsr,ns,tot,xo(),xn(),mloddsrse) ' ' print out results for each 2x2 table ' ? ?#2, ?" TABLE OR Woolf 95% CL Cornfield 95% CL B-D X";chr$(253) ?#2," TABLE OR Woolf 95% CL Cornfield 95% CL B-D X";chr$(253) ?" ----- -- ---------------- ------------------ -------" ?#2," ----- -- ---------------- ------------------ -------" for i=1 to ns oddsr=(xo(1,1,i)+0.5)*(xo(2,2,i)+0.5)/(xo(1,2,i)+0.5)/(xo(2,1,i)+0.5) oddse=sqr(1/(xo(1,1,i)+0.5)+1/(xo(2,2,i)+0.5)+_ 1/(xo(1,2,i)+0.5)+1/(xo(2,1,i)+0.5)) ll=exp(log(oddsr)-1.96*oddse) ul=exp(log(oddsr)+1.96*oddse) if ul>999.99 then ul=999.99 ?" # ";:?using"## ###.##";i;oddsr; ?#2," # ";:?#2,using"## ###.##";i;oddsr; ?using " ###.## to ###.##";ll;ul; ?#2,using " ###.## to ###.##";ll;ul; call corny(xo(1,1,i),xo(1,2,i),xo(2,1,i),xo(2,2,i),ul,ll,1) ?using " ###.## to ###.##";ll;ul; ?#2,using " ###.## to ###.##";ll;ul; call fit(xo(1,1,i),xo(1,2,i),xo(2,1,i),xo(2,2,i),mhoddsr,chitab) ?using " ####.#";chitab ?#2,using " ####.#";chitab next i df=ns-1 ' ' print out pooled results ' ? ?#2, ?"Common Odds Ratio (M-H) ="; ?#2,"Common Odds Ratio (M-H) ="; ?using" ###.##";mhoddsr; ?#2,using" ###.##";mhoddsr; lnodd=log(sumr/sums) selmh=sqr(sumpr/2/sumr/sumr+sumpsqr/2/sumr/sums+sumqs/2/sums/sums) ?" [G-R 95% CI=";:?using" ###.##";EXP(lnodd-1.96*selmh); ?#2," [G-R 95% CI=";:?#2,using" ###.##";EXP(lnodd-1.96*selmh); ?", ";:? USING"###.##]";EXP(lnodd+1.96*selmh) ?#2,", ";:?#2, USING"###.##]";EXP(lnodd+1.96*selmh) ?"Common Odds Ratio (MLE) ="; ?#2,"Common Odds Ratio (MLE) ="; ?using" ###.##";mloddsr; ?#2,using" ###.##";mloddsr; lnodd=log(mloddsr) mloddsrse=log(mloddsrse) ?" [Jack 95% CI=";:?using" ###.##";EXP(lnodd-1.96*selmh); ?#2," [Jack 95% CI=";:?#2,using" ###.##";EXP(lnodd-1.96*selmh); ?", ";:? USING"###.##]";EXP(lnodd+1.96*selmh) ?#2,", ";:?#2, USING"###.##]";EXP(lnodd+1.96*selmh) ? ?#2, ?"LR X";chr$(253);" for homog. of ORs =";:?using"####.##";chisq; ?#2,"LR X";chr$(253);" for homog. of ORs =";:?#2,using"####.##";chisq; ?using" (df=##)";df; ?#2,using" (df=##)";df; call chiprob(chisq,df) chi3=chisq ' ' B-D homog chi-square ' chisq=0 chitab=0 for i=1 to ns call fit(xo(1,1,i),xo(1,2,i),xo(2,1,i),xo(2,2,i),mhoddsr,chitab) chisq=chisq+chitab next i ?"BD X";chr$(253);" for homog. of ORs =";:?using"####.##";chisq; ?#2,"BD X";chr$(253);" for homog. of ORs =";:?#2,using"####.##";chisq; ?using" (df=##)";df; ?#2,using" (df=##)";df; call chiprob(chisq,df) chisq=0 for i=1 to ns: for j=1 to 2: for k=1 to 2 ttot=xo(1,1,i)+xo(1,2,i)+xo(2,1,i)+xo(2,2,i) if xo(j,k,i)<>0 then chisq=chisq+2*xo(j,k,i)*log(xo(j,k,i)_ *ttot/(xo(j,1,i)+xo(j,2,i))/(xo(1,k,i)+xo(2,k,i))) next: next: next chisq=chisq-chi3 ?"Cond. LR X";chr$(253);" (comm OR=1) =";:?using"####.##";chisq; ?#2,"Cond. LR X";chr$(253);" (comm OR=1) =";:?#2,using"####.##";chisq; ?" (df= 1)"; ?#2," (df= 1)"; df=1 call chiprob(chisq,df) chisq=chimh*chimh/var ?"Mantel-H X";chr$(253);" (comm OR=1) =";:?using"####.##";chisq; ?#2,"Mantel-H X";chr$(253);" (comm OR=1) =";:?#2,using"####.##";chisq; ?" (df= 1)"; ?#2," (df= 1)"; df=1 call chiprob(chisq,df) return ' 21000 ' calculator def fngetnum%(a$) if poi%<81 and ((a$=>"0" and a$=<"9") or a$=".") then fngetnum%=1 else fngetnum%=0 end if end def def fnleftb%(s) shared flag%(),stack() local j,t t=0 for j=9 to 22 if stack(s)=j and flag%(s)=1 then t=1 next j fnleftb%=t end def def fnritb%(s) shared flag%(),stack() if stack(s)=8 and flag%(s)=1 then fnritb%=1 else fnritb%=0 end if end def sub eval(oppos) shared flag%(),stack() ' ' will evaluate expression and replace with answer, pulling up stack ' if stack(oppos)=2 then answer=stack(oppos-1)+stack(oppos+1) call rep(oppos,2,answer) elseif stack(oppos)=3 and flag%(oppos-1)=2 then answer=stack(oppos-1)-stack(oppos+1) call rep(oppos,2,answer) elseif stack(oppos)=4 then answer=stack(oppos-1)*stack(oppos+1) call rep(oppos,2,answer) elseif stack(oppos)=5 then if stack(oppos+1)<>0 then answer=stack(oppos-1)/stack(oppos+1) else ?"Division by zero" ?#2,"Division by zero" answer=1e99 end if call rep(oppos,2,answer) elseif stack(oppos)=7 then answer=stack(oppos-1)^stack(oppos+1) call rep(oppos,2,answer) elseif stack(oppos)=6 then answer=stack(oppos-1)^stack(oppos+1) call rep(oppos,1,answer) end if end sub sub rep(oppos,expsize,answer) ' ' actually puts answer in stack and removes expression ' shared flag%(),stack(),size,pos2% flag%(oppos-1)=2 stack(oppos-1)=answer for i=oppos to size-expsize flag%(i)=flag%(i+expsize) stack(i)=stack(i+expsize) next i size=size-expsize pos2%=pos2%-expsize end sub sub dobrack(pos1%,eflag%) shared flag%(),stack() eflag=0 if stack(pos1%)=9 and flag%(pos1%+1)=2 then answer=stack(pos1%+1) call rep(pos1%+1,2,answer) elseif stack(pos1%)=10 and stack(pos1%+1)>0 then answer=log(stack(pos1%+1)) call rep(pos1%+1,2,answer) elseif stack(pos1%)=11 then answer=exp(stack(pos1%+1)) call rep(pos1%+1,2,answer) elseif stack(pos1%)=12 then answer=sin(stack(pos1%+1)) call rep(pos1%+1,2,answer) elseif stack(pos1%)=13 then answer=cos(stack(pos1%+1)) call rep(pos1%+1,2,answer) elseif stack(pos1%)=14 then answer=tan(stack(pos1%+1)) call rep(pos1%+1,2,answer) elseif stack(pos1%)=15 then answer=fix(stack(pos1%+1)) call rep(pos1%+1,2,answer) elseif stack(pos1%)=16 and stack(pos1%+1)>0 then answer=sqr(stack(pos1%+1)) call rep(pos1%+1,2,answer) elseif stack(pos1%)=17 and abs(stack(pos1%+1)<1) then answer=atn(sqr(stack(pos1%+1)^2/(1-stack(pos1%+1)^2))) call rep(pos1%+1,2,answer) elseif stack(pos1%)=17 and abs(stack(pos1%+1)=1) then answer=sgn(stack(pos1%+1))*1.5707966327 call rep(pos1%+1,2,answer) elseif stack(pos1%)=18 and abs(stack(pos1%+1)<1) then answer=atn(sqr((1-stack(pos1%+1)^2)/stack(pos1%+1)^2)) call rep(pos1%+1,2,answer) elseif stack(pos1%)=18 and abs(stack(pos1%+1)=1) then answer=0 call rep(pos1%+1,2,answer) elseif stack(pos1%)=19 then answer=atn(stack(pos1%+1)) call rep(pos1%+1,2,answer) elseif stack(pos1%)=20 and stack(pos1%+1)<1 then answer=0.5*log((1+abs(stack(pos1%+1)))/(1-abs(stack(pos1%+1)))) call rep(pos1%+1,2,answer) elseif stack(pos1%)=21 then answer=1 for c%=1 to stack(pos1%+1) answer=answer*c% next c% call rep(pos1%+1,2,answer) elseif stack(pos1%)=22 and stack(pos1%+1)<10 then answer=FNzp(stack(pos1%+1)) call rep(pos1%+1,2,answer) else eflag%=1 end if end sub op$(1)="=" op$(2)="+" op$(3)="-" op$(4)="*" op$(5)="/" op$(6)="%" op$(7)="^" op$(8)=")" op$(9)="(" op$(10)="log(" op$(11)="exp(" op$(12)="sin(" op$(13)="cos(" op$(14)="tan(" op$(15)="fix(" op$(16)="sqr(" op$(17)="asin(" op$(18)="acos(" op$(19)="atan(" op$(20)="inht(" op$(21)="fact(" op$(22)="zpro(" eflag%=0 randomize timer ' ' main loop ' top: eflag%=0 line input "calc>";lin$ ?#2,"calc> ";lin$ if left$(lin$,1)="q" then return elseif left$(lin$,1)="?" or left$(lin$,1)="h" then ?" --------Calculator---------" ?" = + - * / % ^ ( ) " ?" pi sqr() log() exp() fix() " ?" sin() cos() tan() atan() " ?" asin() acos() inht() rand " ?" fact() zpro() [eval z dist]" ?" assignable variables m1-m9 " ?" format w.d [w=width,d=dec] " ?" quit or q to quit " elseif left$(lin$,6)="format" then poi%=7 w$="":d$="" while mid$(lin$,poi%,1)<>"." and poi%<20 w$=w$+mid$(lin$,poi%,1) poi%=poi%+1 wend poi%=poi%+1 while mid$(lin$,poi%,1)<>" " and poi%<20 d$=d$+mid$(lin$,poi%,1) poi%=poi%+1 wend w=val(w$):d=val(d$) if w>0 and w-d>=0 then format$=string$(w,"#") mid$(format$,w-d,1)="." else format$=string$(w,"#") mid$(format$,1,1)="." end if goto top else if mid$(lin$,1,1)="m" and mid$(lin$,3,1)="=" then poi%=4 assign=val(mid$(lin$,2,1)) else poi%=1 assign=0 end if stapoi%=1 flag%(1)=1 stack(1)=9 maxpoi%=len(lin$)+1 lb%=0 rb%=0 for i=1 to maxpoi% if mid$(lin$,i,1)=")" then incr rb% if mid$(lin$,i,1)="(" then incr lb% next i if rb%<>lb% then ? "Unbalanced parentheses" ?#2,"Unbalanced parentheses" goto top end if while poi%8 and flag%(stapoi%+1)=1 _ and stack(stapoi%+1)=3 then call rep(stapoi%+2,1,-1*stack(stapoi%+2)) if flag%(stapoi%)=2 and flag%(stapoi%+1)=1 _ and stack(stapoi%+1)=6 then call rep(stapoi%+1,1,stack(stapoi%)/100.0) stapoi%=stapoi%+1 wend ' ' get innermost bracket ' while size>1 pos1%=0: pos2%=0: stapoi%=1 while stapoi%=size poi2%=poi2%+1 loop if fnritb%(poi2%)=1 then pos2%=poi2% end if stapoi%=stapoi%+1 wend ' ' evaluate contents ' first order operators by precedence ' while pos2%-pos1%>3 prec=0: oppos=0 for i=pos1%+1 to pos2%-1 if flag%(i)=1 and stack(i)>prec then prec=stack(i):oppos=i next i call eval(oppos) wend call dobrack(pos1%,eflag%) if eflag%=1 then ?"Error in expression" ?#2,"Error in expression" goto top end if wend if assign=0 then if format$="" then ?" =";stack(1) ?#2," =";stack(1) else ?" ="; ?#2," ="; ?using format$;stack(1) ?#2,using format$;stack(1) end if else ?" m"+chr$(48+assign); ?#2," m"+chr$(48+assign); ?"=";stack(1) ?#2,"=";stack(1) vars(assign)=stack(1) end if end if ' ' end main if statement ' goto top ' ' Matched case-control OR and exact 95% CI ' 22000 ' ?"----------------------------------------------------" ?#2,"----------------------------------------------------" ?" Matched case-control OR and Clopper-Pearson 95% CI " ?#2," Matched case-control OR and Clopper-Pearson 95% CI " ?"----------------------------------------------------" ?#2,"----------------------------------------------------" ? ?#2, ?"Enter number of pairs Exposure POS / NEG : ";:INPUT, b ?#2,"Enter number of pairs Exposure POS / NEG : "; b ?"Enter number of pairs Exposure NEG / POS : ";:INPUT, c ?#2,"Enter number of pairs Exposure NEG / POS : "; c call ci95(b,b+c,upp95,low95) ? ?#2, ?using "Odds Ratio = ###.##";b/c ?#2,using "Odds Ratio = ###.##";b/c low95=low95/(1.0d0-low95) if upp95<(b+c) then upp95=upp95/(1.0d0-upp95) if upp95>999.99 then upp95=999.99 else upp95=999.99 end if ?using "Exact 95%CL= ###.## to ###.##";low95;upp95 ?#2,using "Exact 95%CL= ###.## to ###.##";low95;upp95 return ' ' Evaluate z-distribution ' sub zprob(z) local p p=FNzp(z) ?"; tail P=";:? USING "#.###"; p ?#2,"; tail P=";:?#2, USING "#.###"; p end sub ' ' Evaluate z-distribution ' def FNzp(z) local zz#,r#,t#,y#,h#,pi# pi#=3.141592653590 y# =0.471404520791 ' sqr(2)/3 zz#=ABS(z) if z<7 then t#=0 for r#=0 to 12 h#=r#+0.5 t#=t#+exp(-h#*h#/9)*sin(zz#*h#*y#)/h# next r# FNzp=1-2*t#/pi# else FNzp=0 end if end def ' ' Evaluate central chi-square ' ' sub chiprob(chisq,df) local r#,degf#,i#,j#,k#,l#,m#,prob# if chisq=0 then prob#=1:goto fin if chisq>200 then prob#=sqr(2*chisq)-sqr(2*df-1) prob#=0.5+sgn(prob#)*(0.5-FNzp(prob#)/2) ?"; tail P=";:? USING "#.###";prob# ?#2,"; tail P=";:?#2, USING "#.###";prob# exit sub end if r#=1 degf#=df FOR i#=degf# TO 2 STEP -2 r#=r#*i# NEXT i# k#=chisq^(INT((degf#+1)/2))*EXP(-chisq/2)/r# IF INT(degf#/2)=degf#/2 THEN j#=1 else j#=SQR(2/chisq/3.14159265) end if l#=1 m#=1 ' ' main loop ' start: degf#=degf#+2 m#=m#*chisq/degf# IF m#<0.0000001 THEN fin l#=l#+m# GOTO start fin: prob#=1-j#*k#*l#: if prob#<0 then prob#=0 ? "; tail P=";:? USING "#.###";prob# ?#2, "; tail P=";:?#2,using "#.###";prob# end sub ' ' Entry routine for 2x2 tables ' sub twobytwo(a,b,c,d) ? ?#2, ?" D+ D- ?" ------------------- ?" E+ a b ?" E- c d ? ?tab(6);:color 0,7:? "a";:color 15,0:?" "; input;,a ?tab(16);:color 0,7:? "b";:color 15,0:?" "; input,b ?tab(6);:color 0,7:? "c";:color 15,0:?" "; input;,c ?tab(16);:color 0,7:? "d";:color 15,0:?" "; input,d ? ?#2, " Table for analysis " ?#2, " D+ D- " ?#2, " ------------------- " ?#2, " E+";:?#2,using" ##### #####";a;b ?#2, " E-";:?#2,using" ##### #####";c;d ?#2, end sub ' ' Cornfield limits for odds ratio ' based on Ipsen Scand J Soc Med 1984;12:121 ' type=1 odds ratio CL ' type=2 relative risk CL ' type=3 risk difference CL ' sub corny(a,b,c,d,ul,ll,type) local j%,dz#,q#,w#,z#,r#() dim r#(2) if b=0 or c=0 then ul=999.9 exit sub end if z#=0.5 for j%=1 to 2 dz#=1 while abs(dz#)>0.00000001 q#=0 while a+z#<0: z#=z#+a:wend while c-z#<0: z#=z#-c:wend while b-z#<0: z#=z#-b:wend while d+z#<0: z#=z#+d:wend if (a<>0) then q#=q#-a*log(1+z#/a) q#=q#-b*log(1-z#/b)-c*log(1-z#/c) if d#<>0 then q#=q#-d#*log(1+z#/d) w#=a/(a+z#)-b/(b-z#)-c/(c-z#)+d/(d+z#) dz#=(2*q#-3.84146)/2/w# z#=z#+dz# wend if type=1 then r#(j%)=(a+z#)*(d+z#)/(b-z#)/(c-z#) if type=2 then r#(j%)=(a+z#)*(c+d)/(a+b)/(c-z#) if type=3 then r#(j%)=(a+z#)/(a+b)-(c-z#)/(c+d) if type=4 then r#(j%)=((a+z#)*(d+z#)-(b-z#)*(c-z#))/_ sqr((a+b)*(a+c)*(b+d)*(c+d)) z#=-0.5 if a*d/b/c=0 then cfin next j% cfin: ul=r#(1):ll=r#(2) end sub ' Sub Quasisym(n,num(2)) ' ' Caussinus i.p.f algorithm for quasi-symmetry from Haberman(1979) ' local row(),col(),ex0(),ex1(),ex2(),exsum() shared MAXCAT dim row(MAXCAT),col(MAXCAT),ex0(MAXCAT,MAXCAT),_ ex1(MAXCAT,MAXCAT),ex2(MAXCAT,MAXCAT),exsum(MAXCAT) ' ' starting values for iteration and marginal subtotals ' for r=1 to n:for c=1 to n row(r)=row(r)+num(r,c) col(r)=col(r)+num(c,r) ex0(r,c)=1 next c ex0(r,r)=0 next r ' ' main loop ' while abs(ex0(2,1)-ex2(2,1))>0.00001 for r=1 to n:exsum(r)=0:next for r=1 to n:for c=1 to n exsum(r)=exsum(r)+ex0(r,c) next:next for r=1 to n:for c=1 to n if r<>c then ex1(r,c)=ex0(r,c)*(row(r)-num(r,r))/exsum(r) end if next:next for c=1 to n:exsum(c)=0:next for r=1 to n:for c=1 to n exsum(r)=exsum(r)+ex1(c,r) next:next for r=1 to n:for c=1 to n if r<>c then ex2(r,c)=ex1(r,c)*(col(c)-num(c,c))/(exsum(c)) end if next:next for r=1 to n:for c=1 to n if r<>c then ex0(r,c)=ex2(r,c)*(num(r,c)+num(c,r))/(ex2(r,c)+ex2(c,r)) end if next:next wend ' chisq=0:gsq=0:ftsq=0:lowtri=0:upptri=0 for r=1 to n: for c=1 to n if r<>c then ft(c)=sqr(num(r,c))+sqr(num(r,c)+1)-sqr(4*ex0(r,c)+1) chisq=chisq+(num(r,c)-ex0(r,c))^2/ex0(r,c) if num(r,c)>0 then gsq=gsq+2*num(r,c)*log(num(r,c)/ex0(r,c)) ftsq=ftsq+ft(c)*ft(c) re(c)=ex0(r,c) end if next:next df=(n-2)*(n-1)/2 ? ?#2, ? "QuasiSym G";chr$(253);"= ";:?using "###.##";gsq;:chisq=gsq ?#2, "QuasiSym G";chr$(253);"= ";:?#2,using "###.##";gsq;:chisq=gsq call chiprob(chisq,df) end sub ' ' homog chisquare for 2x2 table for given common OR ' sub fit(a,b,c,d,oddsr,chitab) local terma#,termb#,termc#,af#,bf#,cf#,df# if oddsr=1 then terma#=(a*d-b*c)/(a+b+c+d) af#=a-terma# bf#=b+terma# cf#=c+terma# df#=d-terma# else terma#=oddsr-1 termb#=oddsr*(a+a+b+c)+d-a termc#=oddsr*(a+c)*(a+b) termc#=sqr(termb#*termb#-4*terma#*termc#) af#=(termc#+termb#)/2/terma# if af#<=0 or af#>(a+b) or af#>(a+c) then af#=(termb#-termc#)/2/terma# end if bf#=a+b-af# cf#=a+c-af# df#=d-a+af# end if chitab=(a-af#)*(a-af#)*(1/af#+1/bf#+1/cf#+1/df#) end sub ' ' ML fitting no three way interaction model 2x2xK tables ' Sub Itprop(mloddsr,ns,tot,xo(3),xn(3),chisq) shared MAXTAB local it,g,gmean,mtot,onh,lx(),gamma(),h() dim gamma(MAXTAB),h(MAXTAB),lx(2,2,MAXTAB) it=0:g=0:chisq=0 while abs(g)<0.999999 gmean=0 mtot=0 onh=0 it=it+1 for i=1 to ns: for j=1 to 2: for k=1 to 2 if it=1 then lx(j,k,i)=log(xn(j,k,i)) else lx(j,k,i)=log(xn(j,k,i))+(xo(j,k,i)-xn(j,k,i))/xn(j,k,i) end if next: next: next for i=1 to ns gamma(i)=lx(1,1,i)-lx(1,2,i)-lx(2,1,i)+lx(2,2,i) h(i)=1/xn(1,1,i)+1/xn(1,2,i)+1/xn(2,1,i)+1/xn(2,2,i) gmean=gmean+gamma(i)/h(i) onh=onh+1/h(i) next i gmean=gmean/onh for i=1 to ns: for j=1 to 2: for k=1 to 2 xn(j,k,i)=lx(j,k,i)-((-1)^(j+k)*(gamma(i)-gmean))/h(i)/xn(j,k,i) mtot=mtot+exp(xn(j,k,i)) next:next:next g=tot/mtot for i=1 to ns: for j=1 to 2: for k=1 to 2 xn(j,k,i)=g*exp(xn(j,k,i)) next: next: next wend ' ' calculate G**2 and common OR from MLE fitting ' for i=1 to ns: for j=1 to 2: for k=1 to 2 if xo(j,k,i)<>0 then chisq=chisq+2*xo(j,k,i)*log(xo(j,k,i)/xn(j,k,i)) next k: next j: next i ' mloddsr=xn(1,1,1)*xn(2,2,1)/xn(1,2,1)/xn(2,1,1) ' end sub ' ' Jackknife for a standard error for MLE Common OR ' Sub Jack(mloddsr,ns,tot,xo(3),xn(3),mloddsrse) local i,flag,j,k,chisq#,joddsr#(),mj# shared MAXTAB dim joddsr#(2,2,MAXTAB) mloddsrse=0 mj#=0 for i=1 to ns: for j=1 to 2: for k=1 to 2 flag=0 if xo(j,k,i)>0 then xo(j,k,i)=xo(j,k,i)-1:flag=-1 call itprop(joddsr#(j,k,i),ns,tot,xo(),xn(),chisq#) if (flag) then xo(j,k,i)=xo(j,k,i)+1 joddsr#(j,k,i)=tot*mloddsr-(tot-1)*joddsr#(j,k,i) mj#=mj#+xo(j,k,i)*joddsr#(j,k,i) next k: next j: next i mj#=mj#/tot for i=1 to ns: for j=1 to 2: for k=1 to 2 mloddsrse=mloddsrse+xo(j,k,i)*(joddsr(j,k,i)-mj#)*(joddsr(j,k,i)-mj#) next k: next j: next i mloddsrse=sqr(mloddsrse/tot/(tot-1)) mloddsr=mj# end sub