***** 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