Hepatitis
R. K. S. Hankin
To cite the hyper2 package in publications, please use Hankin (2017).
hepatitis_table <- read.table("hepatitis.txt",header=TRUE)
hepatitis_table[,1:3] <- hepatitis_table[,1:3] > 0
wanted <- !(apply(hepatitis_table[,1:3],1,sum) %in% c(0,3))
hepatitis_table <- hepatitis_table[wanted,]
hepatitis_table## P Q E OC LNF
## 2 FALSE FALSE TRUE 63 61
## 3 FALSE TRUE FALSE 55 58
## 4 FALSE TRUE TRUE 18 17
## 5 TRUE FALSE FALSE 69 68
## 6 TRUE FALSE TRUE 17 20
## 7 TRUE TRUE FALSE 21 19W <- hyper2(pnames=c("P","Q","E"))
Wi <- list()
for(i in seq_len(nrow(hepatitis_table))){
jj <- unlist(hepatitis_table[i,1:3])
negative <- names(jj)[jj==0]
positive <- names(jj)[jj==1]
Wi[[i]] <- ggrl(W,positive,negative)
}
hepatitis <- lsl(Wi,powers= hepatitis_table$OC)hepatitis## $suplists
## $suplists[[1]]
## [[1]]
## log( E * (E + P + Q)^-1 * P * (P + Q)^-1)
##
## [[2]]
## log( E * (E + P + Q)^-1 * (P + Q)^-1 * Q)
##
## attr(,"class")
## [1] "list" "suplist"
##
## $suplists[[2]]
## [[1]]
## log((E + P)^-1 * (E + P + Q)^-1 * P * Q)
##
## [[2]]
## log( E * (E + P)^-1 * (E + P + Q)^-1 * Q)
##
## attr(,"class")
## [1] "list" "suplist"
##
## $suplists[[3]]
## [[1]]
## log( E * (E + P)^-1 * (E + P + Q)^-1 * Q)
##
## [[2]]
## log( E * (E + P + Q)^-1 * (P + Q)^-1 * Q)
##
## attr(,"class")
## [1] "list" "suplist"
##
## $suplists[[4]]
## [[1]]
## log((E + P + Q)^-1 * (E + Q)^-1 * P * Q)
##
## [[2]]
## log( E * (E + P + Q)^-1 * (E + Q)^-1 * P)
##
## attr(,"class")
## [1] "list" "suplist"
##
## $suplists[[5]]
## [[1]]
## log( E * (E + P + Q)^-1 * (E + Q)^-1 * P)
##
## [[2]]
## log( E * (E + P + Q)^-1 * P * (P + Q)^-1)
##
## attr(,"class")
## [1] "list" "suplist"
##
## $suplists[[6]]
## [[1]]
## log((E + P + Q)^-1 * (E + Q)^-1 * P * Q)
##
## [[2]]
## log((E + P)^-1 * (E + P + Q)^-1 * P * Q)
##
## attr(,"class")
## [1] "list" "suplist"
##
##
## $powers
## [1] 63 55 18 69 17 21
##
## attr(,"class")
## [1] "lsl"There are several natural BT strengths to consider. First, we could simply count how many times each clinician gives a positive diagnosis:
hepatitis_count <- colSums(sweep(hepatitis_table[,1:3],1,hepatitis_table$OC,"*"))
names(hepatitis_count) <- colnames(hepatitis_table)[1:3]
hepatitis_count <- hepatitis_count/sum(hepatitis_count)
hepatitis_count## P Q E
## 0.35786 0.31438 0.32776Secondly, Zipf:
z <- zipf(3)
names(z) <- c("P","Q","E")
z## P Q E
## 0.54545 0.27273 0.18182Thirdly, equal probabilities:
e <- equalp(W)
e## P Q E
## 0.33333 0.33333 0.33333To find the maximum likelihood estimate we would use:
hepatitis_maxp <- maxp_lsl(hepatitis,startp=indep(hepatitis_count), control=list(trace=110))## Nelder-Mead direct search function minimizer
## function value for initial parameters = 266.484793
## Scaled convergence tolerance is 3.97093e-06
## Stepsize computed as 0.035786
## BUILD 3 267.523381 266.484793
## LO-REDUCTION 5 267.358872 266.484793
## LO-REDUCTION 7 266.834672 266.484793
## HI-REDUCTION 9 266.635763 266.484793
## LO-REDUCTION 11 266.585986 266.484793
## LO-REDUCTION 13 266.502366 266.484793
## HI-REDUCTION 15 266.492180 266.480687
## LO-REDUCTION 17 266.484793 266.480687
## HI-REDUCTION 19 266.480765 266.479273
## HI-REDUCTION 21 266.480687 266.478818
## HI-REDUCTION 23 266.479273 266.478498
## HI-REDUCTION 25 266.478818 266.478494
## REFLECTION 27 266.478498 266.478346
## HI-REDUCTION 29 266.478494 266.478271
## HI-REDUCTION 31 266.478346 266.478247
## HI-REDUCTION 33 266.478271 266.478247
## LO-REDUCTION 35 266.478248 266.478231
## HI-REDUCTION 37 266.478247 266.478225
## HI-REDUCTION 39 266.478231 266.478225
## LO-REDUCTION 41 266.478230 266.478224
## Exiting from Nelder Mead minimizer
## 43 function evaluations used
## Nelder-Mead direct search function minimizer
## function value for initial parameters = 266.478224
## Scaled convergence tolerance is 3.97083e-06
## Stepsize computed as 0.036093
## BUILD 3 267.616723 266.478224
## LO-REDUCTION 5 267.535829 266.478224
## LO-REDUCTION 7 266.919989 266.478224
## LO-REDUCTION 9 266.617501 266.478224
## HI-REDUCTION 11 266.602328 266.478224
## LO-REDUCTION 13 266.535345 266.478224
## LO-REDUCTION 15 266.495953 266.478224
## HI-REDUCTION 17 266.493181 266.478224
## LO-REDUCTION 19 266.485370 266.478224
## LO-REDUCTION 21 266.480583 266.478224
## LO-REDUCTION 23 266.480007 266.478224
## LO-REDUCTION 25 266.479105 266.478224
## LO-REDUCTION 27 266.478482 266.478224
## HI-REDUCTION 29 266.478436 266.478224
## LO-REDUCTION 31 266.478319 266.478224
## LO-REDUCTION 33 266.478253 266.478224
## HI-REDUCTION 35 266.478250 266.478224
## LO-REDUCTION 37 266.478234 266.478224
## Exiting from Nelder Mead minimizer
## 39 function evaluations usedhepatitis_maxp## P Q E
## 0.36093 0.31396 0.32511With this, we are now in a position to compare these four points:
lz <- loglik_lsl(z,hepatitis)
le <- loglik_lsl(e,hepatitis)
lm <- loglik_lsl(hepatitis_maxp,hepatitis)
lc <- loglik_lsl(hepatitis_count,hepatitis)
(results <- c(zipf=lz, equal=le, maxlike=lm, count=lc))## zipf equal maxlike count
## -289.21 -266.96 -266.48 -266.48results <- results - max(results)
results## zipf equal maxlike count
## -22.7281642 -0.4846714 0.0000000 -0.0065688We can use Wilks here to assess the null of hepatitis_count. We
have \(-2\log(\Lambda/\Lambda_0)\sim\chi^2_6\), giving us a \(p\)-value of
pchisq(-2*results[4],df=2,lower.tail=FALSE)## count
## 0.99345Agresti, table 13.1, p542.
Package dataset
Following lines create hepatitis.rda, residing in the data/
directory of the package.
save(hepatitis_table, hepatitis, hepatitis_count, hepatitis_maxp, file="hepatitis.rda")
Hankin, R. K. S. 2017. “Partial Rank Data with the
hyper2 Package: Likelihood Functions for Generalized Bradley-Terry Models.” The R Journal 9 (2): 429–39.