CI.TABLE  Examples.

Author:   Joseph B. Lang,  Department of Statistics and Actuarial Science,  Univ of Iowa, Iowa City, IA 52242 USA

Created:  1/4/09 (Revised 5/18/09)

The examples in this document illustrate how to use the R program  ci.table   (versions 1.1,1.2, and 2.0)  for computing confidence intervals for contingency table parameters.    The examples  herein include those appearing in the first reference below.

Primary References:

Lang, J.B. (2008). "Score and Profile Likelihood Confidence Intervals for Contingency Table Parameters," Statistics in Medicine, 27, 5975-5990.   DOI: 10.1002/sim.3391

Lang, J.B. (2004). "Multinomial Poisson Homogeneous Models for Contingency  Tables," Annals of Statistics, 32, 340-383.

Lang, J.B. (2005). "Homogeneous Linear Predictor Models for Contingency Tables,"  Journal of the American Statistical Association, 100, 121-134.

Availability:  To obtain a free copy of ci.table and accompanying programs, email the author, Joseph B. Lang, at joseph-lang@uiowa.edu.

Disclaimer:  The ci.table code has not been thoroughly tested!

```
###########################################################################
#   BINOMIAL PROBABILITY:  Example 2.1 in Lang, J.B. (2008), Stat Medicine.
###########################################################################
#
#  Observed 1 success in 12 trials.
#  Model: y = (1,11) <- Y ~ mult(12, p1,p2).
#
#  Goal:  Compute an approximate 95% CI for the
#         success probability p1.
#
y <- scan()
1 11

S.fct <- function(p) {
p[1]
}
a1 <- ci.table(y,S.fct,lowerbound=0,upperbound=1,save=T);
a2 <- ci.table(y,S.fct,0,1,type="lr",save=T)
a3 <- ci.table(y,S.fct,0,1,type="pd",save=T)
cbind(a1[c(1:10,20,22)],a2[c(1:10,20,22)],a3[c(1:10,20,22)]) ```
```                    [,1]        [,2]               [,3]
Shat                0.08333333  0.08333333         0.08333333
ase                 0.0797856   0.0797856          0.0797856
lowCI               0.01486509  0.004951072        0.01184658
highCI              0.3538799   0.3185427          0.3437542
cc                  0.95        0.95               0.95
type                "score"     "likelihood ratio" "power divergence"
pdlambda            1           0                  0.6666667
lowWCI              -0.07304355 -0.07304355        -0.07304355
highWCI             0.2397102   0.2397102          0.2397102
calls.to.mph.fit    13          13                 14
time.elapsed        0.06        0.07               0.06```
```> a1\$version
[1] "ci.table (ver 1.1, 2/17/09)"```
```> a4 <- ci.table(y,S.fct,0,1,type="wald",save=T); as.matrix(a4)
Shat                0.08333333
ase                 0.0797856
type                "wald"
lowWCI              -0.07304355
highWCI             0.2397102
version             "ci.table (ver 1.2, 2/27/09)"
time.elapsed        0

#############################################################
#  CONDITIONAL PROBABILITY
#############################################################
#
#  y=(0, 39, 18, 11) <- Y ~ mult(n=68, p1,p2,p3,p4)
#  Goal:  Compute an approximate 95% CI for p1/(p1+p2).
#
y <- scan()
0 39
18 11

S.fct <- function(p) {
p[1]/(p[1]+p[2])
}
ci.table(y,S.fct,lowerbound=0,upperbound=1)```
```                    [,1]
Shat                0
ase                 0.003107673
lowCI               0
highCI              0.08966646
cc                  "95%"
type                "score"
pdlambda            1
lowWCI              -0.006090928
highWCI             0.006090928
lowCI.TestStat      Numeric,0
highCI.TestStat     3.841440
cv                  3.841459
iterlow             0
iterhigh            4
calls.to.mph.fit    9
version             "ci.table (ver 1.1, 2/17/09)"
time.elapsed        0.06            ```
```#############################################################
#  CONDITIONAL PROBABILITY  (Stratified Sampling)
#############################################################
#
#  y = (0, 39 // 18, 11) <- Y ~ prod mult(39, p1, p2 //  29, p3, p4)
#  That is,
#  y <-  Y ~ MP(ν, p| strata=c(1,1,2,2),fixed="all")
#        where   ν = (39, 29).
#
#  Goal:  Compute an approximate 95% CI for p1.
#
y <- scan()
0 39
18 11

strata <- c(1,1,2,2)
S.fct <- function(p) {
p[1]
}
ci.table(y,S.fct,lowerbound=0,upperbound=1,strata=strata)   ```
```                    [,1]
Shat                0
ase                 0.003107673
lowCI               0
highCI              0.08966646
cc                  "95%"
type                "score"
pdlambda            1
lowWCI              -0.006090928
highWCI             0.006090928
lowCI.TestStat      Numeric,0
highCI.TestStat     3.841431
cv                  3.841459
iterlow             0
iterhigh            4
calls.to.mph.fit    9
version             "ci.table (ver 2.0, 5/18/09)"
time.elapsed        0.09

#############################################################
#  DIFFERENCE BETWEEN CONDITIONAL PROBABILITIES
#############################################################
#
#  y=(0, 39, 18, 11) <- Y ~ mult(n=68, p1,p2,p3,p4)
#  Goal:  Compute an approximate 95% CI for p1/(p1+p2) - p3/(p1+p3).
#

y <- scan()
0 39
18 11

S.fct <- function(p) {
p[1]/(p[1]+p[2]) - p[3]/(p[3]+p[4])
}
ci.table(y,S.fct,lowerbound=-1,upperbound=1)```
```                    [,1]
Shat                -0.6206897
ase                 0.09015582
lowCI               -0.7731197
highCI              -0.4400254
cc                  "95%"
type                "score"
pdlambda            1
lowWCI              -0.7973918
highWCI             -0.4439875
lowCI.TestStat      3.841456
highCI.TestStat     3.841459
cv                  3.841459
iterlow             3
iterhigh            3
calls.to.mph.fit    13
version             "ci.table (ver 1.1, 2/17/09)"
time.elapsed        0.17                        ```
```

###################################################################
#  GAMMA VARIANT:  Example 2.3 in Lang, J.B. (2008), Stat Medicine.
###################################################################
#
# y = (25, 25, 12 // 0, 1, 3)  ~ prod Mult(62, p11,p12,p13 // 4, p21,p22,p23)
#
# Goal:  Compute an approximate 95% CI for Gamma* parameter
#        as described in (Lang 2008, Stat Medicine)
#

y <- scan()
25 25 12
0 1 3

strata <- c(1,1,1,2,2,2)

S.fct <- function(p) {
p <- matrix(p, 2,3,byrow=T)
P.case.gt.control <- (p[2,2]+p[2,3])*p[1,1] + p[2,3]*p[1,2]
P.case.lt.control <- p[1,2]*p[2,1] + p[1,3]*(p[2,1]+p[2,2])
P.case.neq.control <- P.case.gt.control + P.case.lt.control
P.case.gt.control/P.case.neq.control
}

ci.table(y,S.fct,0,1,strata=strata)```
```                    [,1]
Shat                0.9358289
ase                 0.06654827
lowCI               0.5250232
highCI              0.98971
cc                  "95%"
type                "score"
pdlambda            1
lowWCI              0.8053967
highWCI             1.066261
lowCI.TestStat      3.841451
highCI.TestStat     3.841457
cv                  3.841459
iterlow             4
iterhigh            3
calls.to.mph.fit    14
version             "ci.table (ver 2.0, 5/18/09)"
time.elapsed        0.37               ```
```
ci.table(y,S.fct,0,1,strata=strata,type="wald")

Shat                0.9358289
ase                 0.06654827
type                "wald"
lowWCI              0.8053967
highWCI             1.066261
version             "ci.table (ver 2.0, 5/18/09)"
time.elapsed        0    ```
```

### Alternative code...

gammastar.fct <- function(p) {
nr <- nrow(p)
nc <- ncol(p)
probC  <- 0;probD <- 0
for (i in 1:(nr-1)) {
for (j in 1:(nc-1)) {
Aij <- 0
for (h in (i+1):nr) {
for (k in (j+1):nc) {
Aij <- Aij + p[h,k]
}
}
probC <- probC + p[i,j]*Aij
}
}
for (i in 1:(nr-1)) {
for (j in 2:nc) {
Aij <- 0
for (h in (i+1):nr) {
for (k in 1:(j-1)) {
Aij <- Aij + p[h,k]
}
}
probD <- probD + p[i,j]*Aij
}
}
probC/(probC+probD)
}

strata <- c(1,1,1,2,2,2)

S.fct <- function(p) {
p <- matrix(p,2,3,byrow=T)
gammastar.fct(p)
}

y <- scan()
25 25 12
0 1 3

a <- ci.table(y,S.fct,0,1,strata=strata,save=T)
b <- ci.table(y,S.fct,0,1,strata=strata,type="pd",save=T)
cbind(a,b)```
```                    a                             b
Shat                0.9358289                     0.9358289
ase                 0.06654827                    0.06654827
lowCI               0.5250232                     0.5367156
highCI              0.98971                       0.9917275
cc                  0.95                          0.95
type                "score"                       "power divergence"
pdlambda            1                             0.6666667
lowWCI              0.8053967                     0.8053967
highWCI             1.066261                      1.066261
lowCI.TestStat      3.841451                      3.841464
highCI.TestStat     3.841457                      3.841417
cv                  3.841459                      3.841459
iterlow             4                             5
iterhigh            3                             3
gSvalues            Numeric,13                    Numeric,14
statvalues          Numeric,13                    Numeric,14
normvalues          Numeric,13                    Numeric,14
warnvalues          Character,13                  Character,14
calls.to.mph.fit    14                            15
version             "ci.table (ver 2.0, 5/18/09)" "ci.table (ver 2.0, 5/18/09)"
time.elapsed        0.48                          1.24
```
```

###################################################################
#  GLOBAL ODDS RATIO
###################################################################
#
# y = (25, 25, 12 // 0, 1, 3)  ~ prod Mult(62, p11,p12,p13 // 4, p21,p22,p23)
#
# Goal:  Compute an approximate 95% CI for a Global Odds Ratio parameter
#        as described in (Lang 2008, Stat Medicine)
#

y <- scan()
25 25 12
0 1 3

strata <- c(1,1,1,2,2,2)

S.fct <- function(p) {
p <- matrix(p, 2,3,byrow=T)
p[1,1]*(p[2,2]+p[2,3])/(p[2,1]*(p[1,2]+p[1,3]))
}
ci.table(y,S.fct,0,NA,strata=strata)```
```Shat                Inf
ase                 1452.629
lowCI               0.6611751
highCI              Inf
cc                  "95%"
type                "score"
pdlambda            1
lowWCI              Inf
highWCI             Inf
lowCI.TestStat      3.84153
highCI.TestStat     Numeric,0
cv                  3.841459
iterlow             4
iterhigh            0
calls.to.mph.fit    8
version             "ci.table (ver 2.0, 5/18/09)"
time.elapsed        0.13

#######################################################
#  DIFFERENCE BETWEEN PRODUCT-MULTINOMIAL PROBABILITIES:
#   Example 2.2 in Lang, J.B. (2008), Stat Medicine.
#######################################################
#
#Source (secondary):  Agresti 2002:65
#Early study of the death penalty in Florida  (Radelet)
#Victim Black...
#  White Defendant  0/9   received Death Penalty
#  Black Defendant  6/103 received Death Penalty
#
#  Model:  y=(0,9 // 6, 97) <- Y ~ prod mult(9, p1, p2 // 103, p3, p4)
#  Goal:  Compute an approximate 95% CI for p1 - p3.
#
y <- c(0,9,6,97)
strata <- c("white","white","black","black")

S.fct <- function(p) {
p[1]-p[3]
}
a1 <- ci.table(y,S.fct,-1,1,strata=strata,save=T)

#Consider the alternative full-multinomial data model...
#   y = (0, 9, 6, 97) <- Y ~ mult(112, p1,p2,p3,p4)
#   Goal:  Compute an approximate 95% CI for p1/(p1+p2) - p3/(p3+p4).
#
#
S.fct <- function(p) {
p[1]/(p[1]+p[2]) - p[3]/(p[3]+p[4])
}

a2 <- ci.table(y,S.fct,-1,1,save=T)

cbind(a1,a2)  ```
```                    a1                            a2
Shat                -0.05825243                   -0.05825243
ase                 0.02534787                    0.02534787
lowCI               -0.1213022                    -0.1213022
highCI              0.2427237                     0.2427237
cc                  0.95                          0.95
type                "score"                       "score"
pdlambda            1                             1
lowWCI              -0.1079333                    -0.1079333
highWCI             -0.008571523                  -0.008571523
lowCI.TestStat      3.841458                      3.841458
highCI.TestStat     3.841449                      3.841449
cv                  3.841459                      3.841459
iterlow             4                             4
iterhigh            5                             5
gSvalues            Numeric,15                    Numeric,15
statvalues          Numeric,15                    Numeric,15
normvalues          Numeric,15                    Numeric,15
warnvalues          Character,15                  Character,15
calls.to.mph.fit    16                            16
version             "ci.table (ver 2.0, 5/18/09)" "ci.table (ver 2.0, 5/18/09)"
time.elapsed        0.29                          0.25

ci.table(y,S.fct,-1,1,type="pd")```
```Shat                -0.05825243
ase                 0.02534787
lowCI               -0.1189673
highCI              0.2169931
cc                  "95%"
type                "power divergence"
pdlambda            0.6666667
lowWCI              -0.1079333
highWCI             -0.008571523
lowCI.TestStat      3.841459
highCI.TestStat     3.841456
cv                  3.841459
iterlow             4
iterhigh            5
calls.to.mph.fit    16
version             "ci.table (ver 2.0, 5/18/09)"
time.elapsed        0.27

##  Alternative (artificial) data that is even more sparse...

ci.table(y=c(0,9,0,97),S.fct,-1,1)```
```                    [,1]
Shat                0
ase                 0.01082153
lowCI               -0.03809444
highCI              0.2991438
cc                  "95%"
type                "score"
pdlambda            1
lowWCI              -0.02120981
highWCI             0.02120981
lowCI.TestStat      3.8415
highCI.TestStat     3.841436
cv                  3.841459
iterlow             3
iterhigh            5
calls.to.mph.fit    15
version             "ci.table (ver 1.1, 2/17/09)"
time.elapsed        0.29

#######################################################################
#  KAPPA COEFFICIENT:  Example 2.4 in Lang, J.B. (2008), Stat Medicine.
#######################################################################
#
#  y = (4,0,0,0,1,0,0,0,15) <-  Y ~ mult(20, p11,p12,...,p33)
#  Goal:  Compute an approximate 95% CI for the
#         unweighted kappa coefficient.
#
y <- c(4,0,0,0,1,0,0,0,15)
S.fct <- function(p) {
p <-  matrix(p,3,3,byrow=T)
s1 <- p[1,1]+p[2,2]+p[3,3]
prow <- apply(p,1,sum)
pcol <- apply(p,2,sum)
s2 <- prow[1]*pcol[1]+prow[2]*pcol[2]+prow[3]*pcol[3]
(s1-s2)/(1-s2)
}
a1 <- ci.table(y,S.fct,-1,1,save=T);
a2 <- ci.table(y,S.fct,-1,1,type="lr",save=T)
a3 <- ci.table(y,S.fct,-1,1,type="pd",save=T)
cbind(a1[c(1:10,20,22)],a2[c(1:10,20,22)],a3[c(1:10,20,22)])  ```
```                    [,1]       [,2]               [,3]
Shat                1          1                  1
ase                 0.06592083 0.06592083         0.06592083
lowCI               0.635835   0.766156           0.665659
highCI              1          1                  1
cc                  0.95       0.95               0.95
type                "score"    "likelihood ratio" "power divergence"
pdlambda            1          0                  0.6666667
lowWCI              0.8707976  0.8707976          0.8707976
highWCI             1.129202   1.129202           1.129202
calls.to.mph.fit    9          8                  10
time.elapsed        4.65       4.61               5   ```
`a1\$version`
```[1] "ci.table (ver 1.1, 2/17/09)"
#######################################################################```

CI.TABLE INPUT and OUTPUT:    (The text below is found in ci.table.Rcode.txt)

```#
# Implementation Instructions:
#
# 1. Source ci.table.Rcode.txt (this file). This creates the two (2) R functions:
#                 ci.table and PD.fct.
#                 Syntax:    R> source("ci.table.Rcode.txt")
#
# 2. Source mph.Rcode.txt.  This creates the nine (9) R functions:
#                 mph.fit,  mph.summary, num.deriv.fct, create.U, create.Z.ZF,
#                 check.homog,  check.HLP, block.fct,  M.fct
#                 Syntax:    R> source("mph.Rcode.txt")
#
# For examples of how to use ci.table, go to
#                 www.stat.uiowa.edu/~jblang/ci.table.documentation
#
#

ci.table <- function(y,S.fct,lowerbound=NA,upperbound=NA,type="score",pdlambda=2/3, cc=0.95,
strata=rep(1,length(y)),fixed.strata="all",
norm.score.conv=1e-5,diffc=0.0001,step=1,app.low=1,app.up=1,
maxiter=25,max.mph.iter=300,epsm=-1,print.output=TRUE,save=FALSE)
{

#ci.table (version 2.0)
#
#   Computes a score, likelihood ratio,  power-divergence, or Wald
#   confidence interval for a contingency table parameter, S(p),
#   where p is a vector of table probabilities.
#   The contingency table counts y are modeled as a realization of
#   a multinomial-Poisson random vector.
#
#   In symbols,
#
#   ci.table  computes a confidence interval for S(p) based on data
#          y <- Y ~ MP(nu, p | strata=s, fixed.strata=F).
#
#   For details on the MP distribution specification, see the mph.fit documentation.
#
#Author:        Joseph B. Lang
#Created:       June 2008,  last updated: 2/27/09,  5/18/09
#
#Primary References:
#
#      Lang, J.B. (2008).  "Score and Profile Likelihood Confidence Intervals for
#          Contingency Table Parameters,"   Statistics in Medicine, 27, 5975-5990.
#          DOI: 10.1002/sim.3391
#
#      Lang, J.B. (2004).  "Multinomial Poisson Homogeneous Models for Contingency
#          Tables,"   Annals of Statistics, 32, 340-383.
#
#      Lang, J.B. (2005).  "Homogeneous Linear Predictor Models for Contingency Tables,"
#          Journal of the American Statistical Association, 100, 121-134.
#
#Subroutines:   ci.table uses PD.fct and  mph.fit (and its subroutines).
#               See the implementation instructions above.
#
#Input:
#
#  Required:   y =vector of table counts
#              S.fct = function object with single argument p, the table probabilities.
#                      S.fct(p) is the parameter of interest.
#
#              Note:  The table counts are [cross-]classification counts.  Generically,
#              y[1] = number of observations with (C = 1),
#              y[2] = number of observations with (C = 2),
#                  etc,
#                     where  C is the [composite] classification variable.
#
#              The [composite] variable C has [joint] distribution characterized by the
#              probability vector P,  defined as P[i] = P(C=i).
#
#              The table probabilities p are interpreted as p = D^{-1}(ZZ'P)P,
#              where P is the vector of joint probabilities and Z is the strata matrix
#              corresponding to strata=s.   Unless strata=rep(1,length(y)),
#              that is, unless a single sample is taken from a single strata,  the table
#              probabilities are conditional probabilities.
#
#  Optional:   lowerbound= smallest possible value for S(p).
#              upperbound= largest possible value for S(p).
#                      Defaults:  lowerbound=NA means no lower bound
#                                 upperbound=NA means no upper bound
#                      !THE USER IS STRONGLY ENCOURAGED TO ENTER NUMERIC VALUES FOR
#                      LOWERBOUND AND UPPERBOUND! (This can mitigate convergence
#                      problems.)
#              type  = interval type, either "score", "lr",  "pd", or "wald"
#                      Default:  type="score"
#                      If type == "lr" the likelihood-ratio (equiv profile likelihood)
#                          confidence interval is computed.
#                      If type == "pd" the power-divergence statistic with lambda=pdlambda
#                          is used to compute the confidence interval.
#                      If type == "wald" only the wald confidence interval is computed.
#                      If type is not any of the four options, the score interval is computed.
#              pdlambda = the parameter used to define the power-divergence statistic.
#                      This parameter is ignored unless type=="pd".  Default: pdlambda=2/3.
#              cc    = the confidence coefficient.  Default:  cc=0.95.
#
#              strata = vector of the same length as y that gives the stratum
#                       membership identifiers.  Default: strata = rep(1,length(y));
#                       i.e. the default is the single stratum (non-stratified) setting.
#                       Examples, strata=A or strata=c(1,1,1,2,2,2,3,3,3) or
#                       strata=paste(sep="","A=",A,",B=",B)
#              fixed.strata = the object that gives information on which stratum have
#                     fixed sample sizes.  It can equal one of the keywords,
#                     fixed.strata = "all" or   fixed.strata = "none",  or it can
#                     be a vector of stratum membership identifiers, e.g.
#                     fixed.strata=c(1,3) or fixed.strata=c("pop1","pop5").
#                     (default: fixed.strata="all")
#                      purposes of computing standard errors for Wald intervals.
#                      Default:  add.to.zero.for.ase = NA, which implies
#              y.eps = amount to initially add to each count in y, as used in mph.fit
#              norm.diff.conv = convergence criterion used in mph.fit.
#                      Default:  norm.diff.conv = 100 (i.e. this criterion is
#                      effectively turned off)
#              norm.score.conv = convergence criterion used in mph.fit.
#              diffc = convergence criterion based on the difference between
#                      current and previous test stat values.
#              step  = step halving parameter as used in mph.fit
#              app.low, app.up   rate parameters that control how fast the
#                      intermediate critical values approach the ultimate critical
#                      value  (which for cc=0.95 is 3.8415).
#              maxiter = maximum number of iterations for the CI root finding algorithm.
#              max.mph.iter = maximum number of iterations as used in mph.fit
#              epsm    =   parameter that controls how far to move from
#                          the empirical estimate Shat (or gShat) when
#                          deriving initial iterates for the CI root
#                          solving algorithm.  For example,  for the upper root,
#                          the inital iterates are gShat+eps, gShat+2*eps, gShat+3*eps,
#                          where, if epsm >= 0, eps=epsm, and if epsm < 0,
#                          eps = 0.5*ase.   Default:  epsm = -1.
#             print.output = logical variable.  If TRUE then print output, otherwise
#                          suppress printing.
#             save   =     logical variable.  If save=FALSE no output is saved.  If
#                          save=TRUE and "a <- ci.table(...)" is submitted then
#                          output is saved in "a".
# Output
#             Shat             =  S.fct(phat),  phat = sample proportions
#             ase              =  approx std error of S.fct(phat)
#             lowCI, highCI    =  endpoints of confidence interval
#             cc               =  nominal confidence coefficient
#             type             =  interval type, either "score" "likelihood ratio"
#                                                or "power divergence"
#             pdlambda         =  power-divergence lambda value used
#             lowWCI, highWCI  =  endpoints of Wald confidence interval
#             lowCI.TestStat,
#             highCI.TestStat  =  test statistic values when the parameter equals the
#                                 lower and upper endpoint of the confidence interval.
#                                 (should be close to the critical value)
#             cv               =  qchisq(cc,1)
#             iterlow, iterhigh=  number of iterations required to find lower and
#                                 upper endpoints
#             gSvalues, statvalues, normvalues, warnvalues = keep track of the g(S(p)) values,
#                                 TestStat values, norm score values, and warning messages that
#                                 result from mph.fit.
#             calls.to.mph.fit =  total number of times restricted MLe had to
#                                 be computed.
#             version          =  version of ci.table that was used
#             time.elapsed     =  time in seconds that it took to compute the confidence interval.

#
# Note 1: To mitigate boundary problems,  if a lower and/or upper bound are input, then
#         a transformed version of the parameter that can take values in -Inf to +Inf
#         is used.   After  finding confidence bounds for the transformed parameter,
#         the endpoints are transformed back to the original scale.
# Note 2: When Shat = - or +Inf,  setting y.eps equal to some positive constant can
#         lead to convergence of the algorithm.  (To make sure the original counts
#         are ultimately used to find the ci bounds, the program should be run
#         again with a very different y.eps value.)
# Note 3: Fine tune with step, app.up, app.low, y.eps, and epsm values.  If
#         lots of 0 counts setting step to a smaller value and perhaps y.eps equal
#         to some positive number may help with convergence.  If TestStat values are
#         sensitive to small changes setting app.low and app.up to smaller values
#         like 0.5 or 0.25 may help.
#
#  ```

Last Updated:  05/18/2009 02:13:16 PM,  Joseph B. Lang

`                       `
`                   `
```
```