//Rate Ratio and SMR analysis on the Internet.
//copyright: Daan G Uitenbroek PhD

function truncate(real,places)  // truncates 'real' numbers n 'places' behind the decimal dot
{
var multi=Math.pow(10,places);
return (Math.round(real*multi)/multi)
}

function cut(real,lenght)  // changes scientific into normal notation and makes to size
{
var str=""+real;
var r="0.";
var pt=str.indexOf("e");
if (pt>0 ) 
    {
     var a=str.substring(0,1)+str.substring(2,pt+1);
     for (var index=1;index<=(str.substring(pt+2,pt+5)*1-1);index++)
        r+=0;
     var st=r+a;
      for (var index=st.length;index<=lenght;index++)
         st+="0";
     return (st.substring(0,lenght));
    };
for (var index=str.length;index<=lenght;index++)
      str+="0";
pt=str.indexOf(".");
if (pt==0) str="0"+str           //solves that newer versions of Netscape give inconsistent numeric notations
return str.substring(0,lenght);
}

function format(integer,fieldlenght) // places leading spaces
{
var str=""+integer;
var r="";
for (var index=1;index<fieldlenght-str.length+1;index++)
 r+=" ";
str=r+str;
return str;
}

function currentCI()
{
this.interval="95%"
this.intzval=1.96
this.intcval=3.8416
this.intup=0.975
this.intlow=0.025
}

function CI(length)
{
this.length=length
   for (var j=1;j<=length;j++)
      this[j]=new currentCI();

this[1].interval="80%"
this[1].intzval=1.28216
this[1].intcval=1.642
this[1].intup=0.90
this[1].intlow=0.10
this[2].interval="90%"
this[2].intzval=1.64485
this[2].intcval=2.706
this[2].intup=0.95
this[2].intlow=0.05
this[3].interval="95%"
this[3].intzval=1.96
this[3].intcval=3.8416
this[3].intup=0.975
this[3].intlow=0.025
this[4].interval="99%"
this[4].intzval=2.57583
this[4].intcval=6.635
this[4].intup=0.995
this[4].intlow=0.005
}

function changeCI(prop)
{
currentoption=prop.selectedIndex+1;
CIcurrent=allCI[currentoption];
}

//following are statements to innitialize the confidence interval
currentoption=3  //3 is for 95% CI
CIcurrent=new currentCI()
allCI=new CI(4)  //because we allow for four different Confidence Intervals

function pois(u,x,q)
{
var writeln="\n";
if (navigator.appVersion.lastIndexOf('Win') != -1) writeln="\r\n";
p0=Math.exp(-u)

tot=1-p0;
document.form1.outp.value+="             POISSON PROBABILITIES"+writeln;
document.form1.outp.value+="           single		cumulative"+writeln;
if (x-3<1)
   document.form1.outp.value+=
        "p(  0): "+cut(p0,10)+";	      p(>  0): "+cut(tot,10)+writeln;
for (index=1;index<=x;index++) 
    {var indextext=""+index;
     var explorer_correct=indextext.length;   //these last two statements correct for microsofts proportional spacing
     p0*=u/index; 
     tot-=p0; 
     if (index>x-3)
     document.form1.outp.value+= 
        "p("+format(index,3)+"): "+cut(p0,11-explorer_correct)+";	      p(>"+format(index,3)+"): "+cut(tot,10)+writeln; 
    }
document.form1.outp.value+=writeln;

rateratio(x,u,q);
}

function fieller(k1,k21,k22)
{
writeln="\n";
if (navigator.appVersion.lastIndexOf('Win') != -1) writeln="\r\n";

var k13=k1/k21
var k14=k21*(k21-CIcurrent.intcval)
var k15=-2*k1*(k21-(k22*CIcurrent.intcval))
var k16=k1*(k1-CIcurrent.intcval)

var k17=Math.sqrt((k15*k15)-(4*k14*k16))
document.form1.outp.value+=0.5*(-k15-k17)/k14+ " <RR< "+-0.5*(k15-k17)/k14 //+"; Chi2 "+CIcurrent.intcval
}

function c2prob(chi2,df)  
{
function odd(number)  
{  
var isodd=0*1;  
if (Math.round(number/2+0.00001)==number/2)  isodd=1*1;  
return isodd  
}  
  
var add=0;  
var mult=1;  
if (odd(df)<0.5)  
 {  
      add=1;  
      mult=Math.sqrt(2/chi2/Math.PI);  
   }  
var denominator=1;  
for (var index=df;index>1;index-=2)  
         denominator=denominator*index;  
var
numerator=Math.pow(chi2,((df*1+add*1)/2))*Math.exp(chi2*-1/2); 

sum=1*1;  
m=1*1;  
  
for (var count=df*1+2*1;m>0.00000000001&&count<20000;count+=2) 

   {  
      m=(m*chi2)/count;  
      sum=sum*1+m*1;  
   }  
return (1-(sum*mult*numerator/denominator));  
}  
  
function invertChi2(p,df) 
{   
var minchisq=0.0*1 
var maxchisq=2000*1 

var chisqval= (df/2) / Math.sqrt (p); 
for (var index=1;maxchisq-minchisq > 0.0000001;index+=1) 
  { 
   if (c2prob(chisqval,df)<p)  
                 maxchisq=chisqval 
   else 
                  minchisq=chisqval 
   chisqval=(maxchisq + minchisq)*0.5 
   } 
return chisqval
} 


function rateratio(Obs,Exp,q)
{
writeln="\n";
if (navigator.appVersion.lastIndexOf('Win') != -1) writeln="\r\n";

var normdev=Math.abs((Obs-Exp)/Math.sqrt(Exp))

var smr=Obs/Exp
var varsmr=Math.sqrt(1/Obs)

var logrr=Math.log(smr)
var selogrr=Math.sqrt((1/Obs)+(1/Exp))

document.form1.outp.value+="Obs= "+Obs+"; Exp= "+Exp+writeln+"Rate Ratio (RR)= "+truncate(smr,4)+"; Std Dev: "+truncate(varsmr,4)+writeln+writeln

document.form1.outp.value+="Exact "+CIcurrent.interval+" Poisson confidence interval: "+writeln
if (Obs>120) document.form1.outp.value+="(Maybe too many cases for reliable estimate)"+writeln
document.form1.outp.value+=(0.5*invertChi2(CIcurrent.intup,(Obs*2))/Exp)+" <RR< "
                        +(0.5*invertChi2(CIcurrent.intlow,(Obs*2+2))/Exp)+writeln


document.form1.outp.value+=writeln+"Estimates with error in expected: "+writeln

var lowerbino=Ibeta(CIcurrent.intlow,Obs,Exp*1+1)
lowerbino=lowerbino/(1-lowerbino)
var upperbino=Ibeta(CIcurrent.intup,Obs*1+1,Exp)
upperbino=upperbino/(1-upperbino)
document.form1.outp.value+="1) The exact "+CIcurrent.interval+" confidence interval: "+writeln
if (Obs>500) document.form1.outp.value+="(Maybe too many cases for reliable estimate)"+writeln
document.form1.outp.value+=lowerbino+"  <RR<  "+upperbino+writeln


document.form1.outp.value+="2) The Fieller "+CIcurrent.interval+" confidence interval: "
                         +writeln
                          fieller(Obs,Exp,q)

var pplow=Obs*Math.pow((1-(1/(9*Obs))-CIcurrent.intzval/Math.sqrt(9*Obs)),3)/Exp
var pphigh=(Obs*1+1*1)*Math.pow((1-(1/(9*(Obs*1+1*1)))+(CIcurrent.intzval/Math.sqrt(9*(Obs*1+1*1)))),3)/Exp

document.form1.outp.value+=writeln+writeln+"Often used approximations:"+writeln
      +CIcurrent.interval+" Walds C.I. "+truncate(Math.exp(Math.log(smr)-(varsmr*CIcurrent.intzval)),4)
      +" <RR< "+truncate(Math.exp(Math.log(smr)+(varsmr*CIcurrent.intzval)),4)+writeln

      +CIcurrent.interval+" Poisson process: "+truncate(pplow,4)
      +" <RR< "+truncate(pphigh,4)+writeln

      +CIcurrent.interval+" Logarithmic Transformation: "+truncate(Math.exp(logrr-(selogrr*CIcurrent.intzval)),4)
      +" <RR< "+truncate(Math.exp(logrr+(selogrr*CIcurrent.intzval)),4)+writeln


      +CIcurrent.interval+" Square-root Transformation: "
      +truncate((Math.sqrt(Obs)-0.5*CIcurrent.intzval)*(Math.sqrt(Obs)-0.5*CIcurrent.intzval)/Exp,4)
      +" <RR< "
      +truncate((Math.sqrt(Obs)+0.5*CIcurrent.intzval)*(Math.sqrt(Obs)+0.5*CIcurrent.intzval)/Exp,4)
      +writeln

      +CIcurrent.interval+" Normal Deviate: "
      +truncate((Obs-Math.sqrt(Obs)*CIcurrent.intzval)/Exp,4)+" <RR< "
      +truncate((Obs*1+Math.sqrt(Obs)*CIcurrent.intzval)/Exp,4)+writeln

 //     +"Normal Dev (meth Rothman)= "+truncate(normdev,2)
 //     +"; probv: "+truncate(tprob(normdev,4000),5)+writeln
 //     +"log ND (conservative)= "+truncate((logrr/selogrr),4)
 //     +"; pr. LND: "+truncate(tprob((logrr/selogrr),4000),5)+writeln
}

function Ibeta(x,p,q)
{
function invertbino(p,x,n) 
{
function bino2(m,k,nk)
{
var p0=Math.exp(Math.log(1-m)*nk)    //(1/m)^nk
var tot=1-p0;
for (var index=0;index<(k-1);index++) 
     {p0*=((nk-index)/(index*1+1*1))*(m/(1-m));
     tot-=p0 
     }
return tot;
}

var minu=0.0*1 
var maxu=1.0*1 

var mean=x/n;   
for (var index=1;maxu-minu > 0.0000001;index+=1) 
  { 
   if (bino2(mean,x,n)>p)  
                 maxu=mean 
   else 
                  minu=mean 
   mean=(maxu*1 + minu*1)*0.5 
   } 
return mean
} 

return invertbino(x,p,(p*1+q*1-1))
}


function ClearOutp() 
{
var writeln="\n";
if (navigator.appVersion.lastIndexOf('Win') != -1) writeln="\r\n";
document.form1.outp.value ="cleared"+writeln;
}