***** SAS IML codes *******;

proc iml;

start cubic(a, b, c, d);
 l=.; u=.;
  if (a = 0) then do;
         if c**2 - 4*b*d > 0 then do;
           l = (-c - sqrt(c**2 - 4*b*d))/(2*b);
           u = (-c + sqrt(c**2 - 4*b*d))/(2*b);
         end;
  end;
  else do;
   a1 = b/a; a2=c/a; a3 = d/a;
   r = a1*a2/6 - a3/2 - a1**3/27;
   q = a2/3 - a1**2/9;
      if q < 0 then do;
        f= r/sqrt(-q**3);
        if abs(f) <=1 then do;
         *if f >1 then f = 1;
         *if f <-1 then f =-1;
         theta = arcos(f);
         pi = 3.1415926;
         k = 2*sqrt(-q)*cos(theta/3)  -a1/3;
         l = 2*sqrt(-q) * cos ( (theta+2*pi)/3  ) - a1/3;
         u = 2*sqrt(-q) * cos ( (theta+4*pi)/3  ) - a1/3;
         end;
     end;
 end;
return(l||u||k);

finish cubic;



start interval4ICC(alpha, data);
 crit = probit( 1- alpha/2);
 chi = crit**2;
 yi = data[,1];
 ni = data[,2];

 bigN = sum(ni);
 piest= sum(yi)/bigN;
 piHat = piest;
 index = piest * (1 - piest);
 k = nrow(ni); * 26; *armclusters;


****** Fleiss-Cuzick **********;
rhoHat = 1 - sum ( yi#(ni-yi)/ni )/((bigN- k)* piest*(1-piest));


part1 = ((1/index - 6) * sum(1/ni))/(bigN - k)**2 + ((2 *bigN
                 + 4 * k - k/index) * k)/(bigN * (bigN - k)**2);
part2 = sum(ni##2)/(bigN**2 * index)-
               (3 * bigN - 2 * k) * (bigN - 2 * k) * sum(ni##2)/
               (bigN**2* (bigN - k)**2) - (2 * bigN - k)/(bigN - k)**2;
part3 = (4-1/index)*(sum(ni##2) - bigN)/bigN**2;
vFC  =  (1 - rhoHat) * (part1 + part2*rhoHat + part3*rhoHat**2);

** CUBIC ****;
A2 = - part3;
B2 = part3 - part2;
C2 = part2 - part1;
D2 = part1;


A22 =chi*A2;
B22 =chi*B2 - 1;
C22 =chi*C2 + 2* rhoHat;
D22 =chi*D2 - rhoHat**2;
solution2 = cubic(A22, B22, C22, D22);
lower = solution2[,1];
upper = solution2[,2];

  print alpha piHat rhoHat lower upper;

finish interval4ICC;


**** Liang et al (1992) ****;

yi = j(36,1,0)//J(12,1,1)//j(15,1,0)//J(7,1,1)//j(1,1,2)
      //j(5,1,0)//j(7,1,1)//j(3,1,2)//J(2,1,3)
      //J(3,1,0)//J(3,1,1)//J(1,1,2)
      //{0,2,3,4,6};
ni = j(48,1,1)//J(23,1,2)//J(17,1,3)//J(7,1,4)//J(5,1,6);
data = yi||ni;

print 'COPD data (Liang et al, 1992)';

run interval4ICC(0.05, data);

* Fleiss example;
yi = {6 0 0 3 0 0 1 4 5 4 0 5 3 0 1 0 4 0 5 1 4 4 1 0 4 0};
ni = {6 3 5 6 6 4 6 6 6 6 6 6 5 5 4 6 6 3 6 3 6 5 6 4 6 6};
 
yi = t(yi);
ni = t(ni);

print 'Fleiss data (Lipsitz et al, 1994)';
data =yi||ni;
run interval4ICC(0.05, data);

quit;


OUtput:



                            COPD data (Liang et al, 1992)


                      ALPHA     PIHAT    RHOHAT     LOWER     UPPER

                       0.05 0.2955665 0.1800851 0.0079251 0.4019612


                          Fleiss data (Lipsitz et al, 1994)


                      ALPHA     PIHAT    RHOHAT     LOWER     UPPER

                       0.05 0.4014599 0.4094969 0.2102443 0.6209409




##### S-Plus Code




Interval4ICC _ function(alpha=0.05, yi, ni)
{
chi _ qchisq(1-alpha, 1)
pi.Hat _ sum(yi)/sum(ni)
index _ pi.Hat * (1 - pi.Hat)
k _ length(ni)
bigN _ sum(ni)
rho.Hat _ 1 - sum ( yi* (ni-yi)/ni )/((bigN- k)* pi.Hat*(1-pi.Hat));

part1 _ ((1/index - 6) * sum(1/ni))/(bigN - k)^2 + ((2 *bigN
                 + 4 * k - k/index) * k)/(bigN * (bigN - k)^2);
part2 _ sum(ni^2)/(bigN^2 * index)-
               (3 * bigN - 2 * k) * (bigN - 2 * k) * sum(ni^2)/
               (bigN^2* (bigN - k)^2) - (2 * bigN - k)/(bigN - k)^2;
part3 _ (4-1/index)*(sum(ni^2) - bigN)/bigN^2;

vFC  _  (1 - rho.Hat) * (part1 + part2*rho.Hat + part3*rho.Hat^2);

A  _ - part3;
B  _ part3 - part2;
C  _ part2 - part1;
D  _ part1;


A1 _ chi*A;
B1 _chi*B - 1;
C1 _chi*C + 2* rho.Hat;
D1 _chi*D - rho.Hat^2;

root _ polyroot(c(D1, C1, B1, A1))

return(pi.Hat, rho.Hat, root)
 
} 



Example
 
fleiss.y _c(6, 0, 0, 3, 0, 0, 1, 4, 5, 4, 0, 5, 3, 0, 1, 0, 4, 0, 5, 1, 4, 4, 1, 0, 4, 0)
fleiss.n _c(6, 3, 5, 6, 6, 4, 6, 6, 6, 6, 6, 6, 5, 5, 4, 6, 6, 3, 6, 3, 6, 5, 6, 4, 6, 6)


Interval4ICC(yi=fleiss.y,ni=fleiss.n)


Output:

$pi.Hat:
[1] 0.4014599

$rho.Hat:
[1] 0.4094969

$root:
[1]  0.2102432-4.038968e-028i  0.6209433+4.102077e-028i
[3] 58.3419310-6.310887e-030i