R/InterfaceOneSampleIR.R
In bjbiggerstaff/PooledInfRate: Estimation for Pooled or Group Testing
"pIR" <- function(x,...){
UseMethod("pIR")
}
"pIR.default" <- function(x,m,n=rep(1,length(x)), group,
pt.method = c("firth","gart","bc-mle","mle","mir"),
ci.method = c("skew-score","bc-skew-score","score","lrt","wald","mir"),
scale=1, alpha=0.05, tol=.Machine$double.eps^0.5, ...) {
call <- match.call()
call[[1]] <- as.name("pIR")
xmn.nomiss <- (!is.na(x) & !is.na(m) & !is.na(n))
x <- x[xmn.nomiss]
m <- m[xmn.nomiss]
n <- n[xmn.nomiss]
mn.nozero <- (m>0 & n>0)
x <- x[mn.nozero]
m <- m[mn.nozero]
n <- n[mn.nozero]
if(!missing(group)){
group.var <- "Group" #deparse(substitute(group))
xmng.nomiss <- (!is.na(x) & !is.na(m) & !is.na(n) & !is.na(group))
group <- group[xmng.nomiss]
groups <- unique(group)
nGroups <- length(groups)
group.dat <- data.frame(Group=group)
} else {
xmng.nomiss <- (!is.na(x) & !is.na(m) & !is.na(n))
group <- rep("SINGLE",length(x[xmng.nomiss]))
groups <- unique(group)
nGroups <- length(groups)
group.var <- ""
group.dat <- data.frame(Group=rep(1,length(x[xmng.nomiss])))
}
ans <- vector(mode="list",length=nGroups)
for(i in 1:nGroups){
ans[[i]] <- pooledBin.fit(x[group==groups[i]],
m[group==groups[i]],
n[group==groups[i]],
pt.method=pt.method,
ci.method=ci.method,
scale=scale,alpha=alpha,tol=tol)
}
ans.lst <- ans
ans <- cbind(group.dat[!duplicated(group),],
data.frame(
P = rep(0, nGroups),
Lower = rep(0, nGroups),
Upper = rep(0, nGroups),
Scale = rep(1, nGroups))
)
#ans <- data.frame(Group = groups,
# P = rep(0, nGroups),
# Lower = rep(0, nGroups),
# Upper = rep(0, nGroups),
# Scale = rep(1, nGroups))
names(ans)[1] <- ifelse(group.var != "",group.var,"Group")
#if (group.var != "") names(ans)[1] <- group.var
for (i in 1:nGroups) ans[i, 2:5] <- ans.lst[[i]]$scale * c(ans.lst[[i]]$p, ans.lst[[i]]$lcl, ans.lst[[i]]$ucl, 1)
if (all(ans$Scale == 1)) ans$Scale <- NULL
if(nGroups == 1) ans$Group <- NULL
# fix row names after subsetting
group.dat <- as.data.frame(group.dat[!duplicated(group),])
rownames(group.dat) <- 1:nrow(group.dat)
ans <- structure(ans, class = "pooledBin",fullList = ans.lst,
group.names = groups, group.var = group.var,n.groups = nGroups,
group.dat = group.dat,
scale=scale,call=call)
ans
}
# "pIR.formula" <- function(x, data,
# pt.method = c("firth","gart","bc-mle","mle","mir"),
# ci.method = c("skew-score","bc-skew-score","score","lrt","wald","mir"),
# scale=1, alpha=0.05, tol=.Machine$double.eps^0.5, ...){
# call <- match.call()
# call[[1]] <- as.name("pIR")
#
# if(missing(data))
# data <- environment(x)
#
# vars <- pooledBinParseFormula(x, data)
# if(any(sapply(vars[1:3],length)>1)) stop("only variable names permitted in formula; perhaps use the default call")
# # omit records with missing data -- note, if data contains records missing
# # anywhere (even not in X, M, N, Group, they are omitted), so care should be
# # used in subsetting before the call to be assured of the desired analysis
# # restrict to only those variables needed before subsetting to avoid that issue.
# if(!missing(data)){
# # restrict the data to the variables needed before removing records with any NA
# # --this doesn't help when no data are specified, so that's dealt with below
# data <- data[,unlist(vars)]
# data <- na.omit(data)
# }
#
# # retrieve values from data using the character name
# # use the eval(parse(text=XX), data) construct in case data is left unspecified,
# # and because we have the names of the variables available
# # from pooledBinParseFormula
# x <- eval(parse(text=vars$x), data)
# m <- eval(parse(text=vars$m), data)
# if(!is.null(vars$n))
# n <- eval(parse(text=vars$n), data)
# else n <- rep(1,length(x)) # default n
#
# if(!is.null(vars$group[1])){
# # NEW as of 5/24/2022 to allow multiple grouping variables
# n.group.var <- length(vars$group)
# group.dat <- as.data.frame(matrix(nrow=length(x),ncol=n.group.var))
# names(group.dat) <- vars$group
#
# for(i in 1:n.group.var) group.dat[,i] <- eval(parse(text=vars$group[i]),data)
#
# #print(interaction(group.dat,drop=TRUE)) # drop unused levels
# #return(group.dat)
#
# # end of NEW
#
# #group <- eval(parse(text=vars$group), data)
# group <- interaction(group.dat,drop=TRUE)
# groups <- unique(group)
# nGroups <- length(groups)
# group.var <- vars$group
# } else {
# group <- rep("SINGLE",length(x))
# groups <- unique(group)
# nGroups <- 1
# group.var <- ""
# group.dat <- data.frame()
# }
#
# # restrict to data with no missing x, m, n, group
# if(missing(data)){
# mn.nozero <- (m>0 & n>0)
# x <- x[mn.nozero]
# m <- m[mn.nozero]
# n <- n[mn.nozero]
#
# if(nGroups == 1){
# xmng.nomiss <- (!is.na(x) & !is.na(m) & !is.na(n))
# } else {
# xmng.nomiss <- (!is.na(x) & !is.na(m) & !is.na(n) & !is.na(group))
# }
# x <- x[xmng.nomiss]
# m <- m[xmng.nomiss]
# n <- n[xmng.nomiss]
#
# group <- group[xmng.nomiss]
# groups <- unique(group)
# nGroups <- length(groups)
# }
#
# mn.nozero <- (m>0 & n>0)
# x <- x[mn.nozero]
# m <- m[mn.nozero]
# n <- n[mn.nozero]
# group <- group[mn.nozero]
# groups <- unique(group)
# nGroups <- length(groups)
#
# #if(nGroups > 1){
# ans <- vector(mode="list",length=nGroups)
# for(i in 1:nGroups){
# ans[[i]] <- pooledBin.fit(x[group==groups[i]],
# m[group==groups[i]],
# n[group==groups[i]],
# pt.method=pt.method,
# ci.method=ci.method,
# scale=scale,alpha=alpha,tol=tol)
# }
# names(ans) <- groups
# ans.lst <- ans
#
# # ans <- data.frame(Group = groups,
# # P = rep(0, nGroups),
# # Lower = rep(0, nGroups),
# # Upper = rep(0, nGroups),
# # Scale = rep(1, nGroups))
# ans <- cbind(group.dat[!duplicated(group),],
# data.frame(
# P = rep(0, nGroups),
# Lower = rep(0, nGroups),
# Upper = rep(0, nGroups),
# Scale = rep(1, nGroups))
# )
# if (!is.null(vars$group)){
# #names(ans)[1] <- paste0(vars$group,collapse=".") # new as of 5/26/2022 for multiple grouping variables
# names(ans)[1:length(vars$group)] <- vars$group
# }
# for (i in 1:nGroups) ans[i,(ncol(group.dat)-1)+ (2:5)] <- ans.lst[[i]]$scale * c(ans.lst[[i]]$p, ans.lst[[i]]$lcl, ans.lst[[i]]$ucl, 1)
# if (all(ans$Scale == 1)) ans$Scale <- NULL
#
# if(nGroups == 1) ans$Group <- NULL
#
# # fix row names after subsetting
# group.dat <- as.data.frame(group.dat[!duplicated(group),])
# rownames(group.dat) <- 1:nrow(group.dat)
#
# ans <- structure(ans, class = "pIR", fullList = ans.lst,
# x.var = vars$x, m.var = vars$m, n.var = vars$n,
# group.names = groups, group.var = group.var,
# group.dat = group.dat,
# n.groups = nGroups, scale=scale,call=call)
# ans
# }
"pIR.formula" <- function(x, data,
pt.method = c("firth","gart","bc-mle","mle","mir"),
ci.method = c("skew-score","bc-skew-score","score","lrt","wald","mir"),
scale=1, alpha=0.05, tol=.Machine$double.eps^0.5, ...){
call <- match.call()
call[[1]] <- as.name("pooledBin")
have.df <- (!missing(data))
if(!have.df){
data <- environment(as.formula(x))
}
#data <- environment(x)
vars <- pooledBinParseFormula(x, data)
if(any(sapply(vars[1:3],length)>1)) stop("only variable names permitted in formula; perhaps use the default call")
# omit records with missing data -- note, if data contains records missing
# anywhere (even not in X, M, N, Group, they are omitted), so care should be
# used in subsetting before the call to be assured of the desired analysis
# restrict to only those variables needed before subsetting to avoid that issue.
#if(!missing(data)){
if(have.df){
# restrict the data to the variables needed before removing records with any NA
# --this doesn't help when no data are specified, so that's dealt with below
data <- data[,unlist(vars)]
data <- na.omit(data)
}
# retrieve values from data using the character name
# use the eval(parse(text=XX), data) construct in case data is left unspecified,
# and because we have the names of the variables available
# from pooledBinParseFormula
x <- eval(parse(text=vars$x), data)
m <- eval(parse(text=vars$m), data)
if(!is.null(vars$n))
n <- eval(parse(text=vars$n), data)
else n <- rep(1,length(x)) # default n
if(!is.null(vars$group[1])){
# NEW as of 5/24/2022 to allow multiple grouping variables
n.group.var <- length(vars$group)
group.dat <- as.data.frame(matrix(nrow=length(x),ncol=n.group.var))
names(group.dat) <- vars$group
for(i in 1:n.group.var) group.dat[,i] <- eval(parse(text=vars$group[i]),data)
#print(interaction(group.dat,drop=TRUE)) # drop unused levels
#return(group.dat)
# end of NEW
#group <- eval(parse(text=vars$group), data)
group <- interaction(group.dat,drop=TRUE)
groups <- unique(group)
nGroups <- length(groups)
group.var <- vars$group
} else {
group <- rep("SINGLE",length(x))
groups <- unique(group)
nGroups <- 1
group.var <- ""
group.dat <- data.frame()
}
# restrict to data with no missing x, m, n, group
#if(missing(data)){
mn.nozero <- (m>0 & n>0)
x <- x[mn.nozero]
m <- m[mn.nozero]
n <- n[mn.nozero]
if(nGroups == 1){
xmng.nomiss <- (!is.na(x) & !is.na(m) & !is.na(n))
} else {
xmng.nomiss <- (!is.na(x) & !is.na(m) & !is.na(n) & !is.na(group))
}
x <- x[xmng.nomiss]
m <- m[xmng.nomiss]
n <- n[xmng.nomiss]
group <- group[xmng.nomiss]
groups <- unique(group)
nGroups <- length(groups)
#}
# mn.nozero <- (m>0 & n>0)
# x <- x[mn.nozero]
# m <- m[mn.nozero]
# n <- n[mn.nozero]
# group <- group[mn.nozero]
# groups <- unique(group)
# nGroups <- length(groups)
#if(nGroups > 1){
ans <- vector(mode="list",length=nGroups)
for(i in 1:nGroups){
ans[[i]] <- pooledBin.fit(x[group==groups[i]],
m[group==groups[i]],
n[group==groups[i]],
pt.method=pt.method,
ci.method=ci.method,
scale=scale,alpha=alpha,tol=tol)
}
names(ans) <- groups
ans.lst <- ans
# ans <- data.frame(Group = groups,
# P = rep(0, nGroups),
# Lower = rep(0, nGroups),
# Upper = rep(0, nGroups),
# Scale = rep(1, nGroups))
ans <- cbind(group.dat[!duplicated(group),],
data.frame(
P = rep(0, nGroups),
Lower = rep(0, nGroups),
Upper = rep(0, nGroups),
Scale = rep(1, nGroups))
)
if (!is.null(vars$group)){
#names(ans)[1] <- paste0(vars$group,collapse=".") # new as of 5/26/2022 for multiple grouping variables
names(ans)[1:length(vars$group)] <- vars$group
}
for (i in 1:nGroups) ans[i,(ncol(group.dat)-1)+ (2:5)] <- ans.lst[[i]]$scale * c(ans.lst[[i]]$p, ans.lst[[i]]$lcl, ans.lst[[i]]$ucl, 1)
if (all(ans$Scale == 1)) ans$Scale <- NULL
if(nGroups == 1) ans$Group <- NULL
# fix row names after subsetting
group.dat <- as.data.frame(group.dat[!duplicated(group),])
rownames(group.dat) <- 1:nrow(group.dat)
ans <- structure(ans, class = "pIR", fullList = ans.lst,
x.var = vars$x, m.var = vars$m, n.var = vars$n,
group.names = groups, group.var = group.var,
group.dat = group.dat,
n.groups = nGroups, scale=scale,call=call)
ans
}
"print.pIR" <- function(x, ...){
print(as.data.frame(unclass(x)),...)
invisible(x)
}
# to build the package, had to use a version of "[.data.frame" without the .Internal(copyDFattr(xx,x)) call
# this just copied x as a list along with the attributes of the data frame
#"[.pIR" <- subsetDF #get("[.data.frame")
#"[.summary.pIR" <- subsetDF #get("[.data.frame")
"[.pIR" <- function (x, i, j, drop = if (missing(i)) TRUE else length(cols) == 1) {
mdrop <- missing(drop)
Narg <- nargs() - !mdrop
has.j <- !missing(j)
if (!all(names(sys.call()) %in% c("", "drop")) &&
!isS4(x))
warning("named arguments other than 'drop' are discouraged")
if (Narg < 3L) {
if (!mdrop)
warning("'drop' argument will be ignored")
if (missing(i))
return(x)
if (is.matrix(i))
return(as.matrix(x)[i])
nm <- names(x)
if (is.null(nm))
nm <- character()
if (!is.character(i) && anyNA(nm)) {
names(nm) <- names(x) <- seq_along(x)
y <- NextMethod("[")
cols <- names(y)
if (anyNA(cols))
stop("undefined columns selected")
cols <- names(y) <- nm[cols]
}
else {
y <- NextMethod("[")
cols <- names(y)
if (!is.null(cols) && anyNA(cols))
stop("undefined columns selected")
}
if (anyDuplicated(cols))
names(y) <- make.unique(cols)
attr(y, "row.names") <- .row_names_info(x, 0L)
attr(y, "class") <- oldClass(x)
return(y)
}
if (missing(i)) {
if (drop && !has.j && length(x) == 1L)
return(.subset2(x, 1L))
nm <- names(x)
if (is.null(nm))
nm <- character()
if (has.j && !is.character(j) && anyNA(nm)) {
names(nm) <- names(x) <- seq_along(x)
y <- .subset(x, j)
cols <- names(y)
if (anyNA(cols))
stop("undefined columns selected")
cols <- names(y) <- nm[cols]
}
else {
y <- if (has.j)
.subset(x, j)
else x
cols <- names(y)
if (anyNA(cols))
stop("undefined columns selected")
}
if (drop && length(y) == 1L)
return(.subset2(y, 1L))
if (anyDuplicated(cols))
names(y) <- make.unique(cols)
nrow <- .row_names_info(x, 2L)
if (drop && !mdrop && nrow == 1L)
return(structure(y, class = NULL, row.names = NULL))
else {
attr(y, "class") <- oldClass(x)
attr(y, "row.names") <- .row_names_info(x,
0L)
return(y)
}
}
xx <- x
cols <- names(xx)
x <- vector("list", length(x))
#x <- .Internal(copyDFattr(xx, x))
x <- as.list(xx)
attributes(x) <- attributes(xx)
oldClass(x) <- attr(x, "row.names") <- NULL
if (has.j) {
nm <- names(x)
if (is.null(nm))
nm <- character()
if (!is.character(j) && anyNA(nm))
names(nm) <- names(x) <- seq_along(x)
x <- x[j]
cols <- names(x)
if (drop && length(x) == 1L) {
if (is.character(i)) {
rows <- attr(xx, "row.names")
i <- pmatch(i, rows, duplicates.ok = TRUE)
}
xj <- .subset2(.subset(xx, j), 1L)
return(if (length(dim(xj)) != 2L) xj[i] else xj[i,
, drop = FALSE])
}
if (anyNA(cols))
stop("undefined columns selected")
if (!is.null(names(nm)))
cols <- names(x) <- nm[cols]
nxx <- structure(seq_along(xx), names = names(xx))
sxx <- match(nxx[j], seq_along(xx))
}
else sxx <- seq_along(x)
rows <- NULL
if (is.character(i)) {
rows <- attr(xx, "row.names")
i <- pmatch(i, rows, duplicates.ok = TRUE)
}
for (j in seq_along(x)) {
xj <- xx[[sxx[j]]]
x[[j]] <- if (length(dim(xj)) != 2L)
xj[i]
else xj[i, , drop = FALSE]
}
if (drop) {
n <- length(x)
if (n == 1L)
return(x[[1L]])
if (n > 1L) {
xj <- x[[1L]]
nrow <- if (length(dim(xj)) == 2L)
dim(xj)[1L]
else length(xj)
drop <- !mdrop && nrow == 1L
}
else drop <- FALSE
}
if (!drop) {
if (is.null(rows))
rows <- attr(xx, "row.names")
rows <- rows[i]
if ((ina <- anyNA(rows)) | (dup <- anyDuplicated(rows))) {
if (!dup && is.character(rows))
dup <- "NA" %in% rows
if (ina)
rows[is.na(rows)] <- "NA"
if (dup)
rows <- make.unique(as.character(rows))
}
if (has.j && anyDuplicated(nm <- names(x)))
names(x) <- make.unique(nm)
if (is.null(rows))
rows <- attr(xx, "row.names")[i]
attr(x, "row.names") <- rows
oldClass(x) <- oldClass(xx)
}
x
}
"[.summary.pIR" <- function (x, i, j, drop = if (missing(i)) TRUE else length(cols) == 1) {
mdrop <- missing(drop)
Narg <- nargs() - !mdrop
has.j <- !missing(j)
if (!all(names(sys.call()) %in% c("", "drop")) &&
!isS4(x))
warning("named arguments other than 'drop' are discouraged")
if (Narg < 3L) {
if (!mdrop)
warning("'drop' argument will be ignored")
if (missing(i))
return(x)
if (is.matrix(i))
return(as.matrix(x)[i])
nm <- names(x)
if (is.null(nm))
nm <- character()
if (!is.character(i) && anyNA(nm)) {
names(nm) <- names(x) <- seq_along(x)
y <- NextMethod("[")
cols <- names(y)
if (anyNA(cols))
stop("undefined columns selected")
cols <- names(y) <- nm[cols]
}
else {
y <- NextMethod("[")
cols <- names(y)
if (!is.null(cols) && anyNA(cols))
stop("undefined columns selected")
}
if (anyDuplicated(cols))
names(y) <- make.unique(cols)
attr(y, "row.names") <- .row_names_info(x, 0L)
attr(y, "class") <- oldClass(x)
return(y)
}
if (missing(i)) {
if (drop && !has.j && length(x) == 1L)
return(.subset2(x, 1L))
nm <- names(x)
if (is.null(nm))
nm <- character()
if (has.j && !is.character(j) && anyNA(nm)) {
names(nm) <- names(x) <- seq_along(x)
y <- .subset(x, j)
cols <- names(y)
if (anyNA(cols))
stop("undefined columns selected")
cols <- names(y) <- nm[cols]
}
else {
y <- if (has.j)
.subset(x, j)
else x
cols <- names(y)
if (anyNA(cols))
stop("undefined columns selected")
}
if (drop && length(y) == 1L)
return(.subset2(y, 1L))
if (anyDuplicated(cols))
names(y) <- make.unique(cols)
nrow <- .row_names_info(x, 2L)
if (drop && !mdrop && nrow == 1L)
return(structure(y, class = NULL, row.names = NULL))
else {
attr(y, "class") <- oldClass(x)
attr(y, "row.names") <- .row_names_info(x,
0L)
return(y)
}
}
xx <- x
cols <- names(xx)
x <- vector("list", length(x))
#x <- .Internal(copyDFattr(xx, x))
x <- as.list(xx)
attributes(x) <- attributes(xx)
oldClass(x) <- attr(x, "row.names") <- NULL
if (has.j) {
nm <- names(x)
if (is.null(nm))
nm <- character()
if (!is.character(j) && anyNA(nm))
names(nm) <- names(x) <- seq_along(x)
x <- x[j]
cols <- names(x)
if (drop && length(x) == 1L) {
if (is.character(i)) {
rows <- attr(xx, "row.names")
i <- pmatch(i, rows, duplicates.ok = TRUE)
}
xj <- .subset2(.subset(xx, j), 1L)
return(if (length(dim(xj)) != 2L) xj[i] else xj[i,
, drop = FALSE])
}
if (anyNA(cols))
stop("undefined columns selected")
if (!is.null(names(nm)))
cols <- names(x) <- nm[cols]
nxx <- structure(seq_along(xx), names = names(xx))
sxx <- match(nxx[j], seq_along(xx))
}
else sxx <- seq_along(x)
rows <- NULL
if (is.character(i)) {
rows <- attr(xx, "row.names")
i <- pmatch(i, rows, duplicates.ok = TRUE)
}
for (j in seq_along(x)) {
xj <- xx[[sxx[j]]]
x[[j]] <- if (length(dim(xj)) != 2L)
xj[i]
else xj[i, , drop = FALSE]
}
if (drop) {
n <- length(x)
if (n == 1L)
return(x[[1L]])
if (n > 1L) {
xj <- x[[1L]]
nrow <- if (length(dim(xj)) == 2L)
dim(xj)[1L]
else length(xj)
drop <- !mdrop && nrow == 1L
}
else drop <- FALSE
}
if (!drop) {
if (is.null(rows))
rows <- attr(xx, "row.names")
rows <- rows[i]
if ((ina <- anyNA(rows)) | (dup <- anyDuplicated(rows))) {
if (!dup && is.character(rows))
dup <- "NA" %in% rows
if (ina)
rows[is.na(rows)] <- "NA"
if (dup)
rows <- make.unique(as.character(rows))
}
if (has.j && anyDuplicated(nm <- names(x)))
names(x) <- make.unique(nm)
if (is.null(rows))
rows <- attr(xx, "row.names")[i]
attr(x, "row.names") <- rows
oldClass(x) <- oldClass(xx)
}
x
}
"summary.pIR" <- function(object, simple = FALSE, ...){
grp.names <- as.character(attr(object,"group.names"))
"sumf" <- function(x) c(x, N = sum(x$n * x$m), NumPools = sum(x$n), NumPosPools = sum(x$x),
PtEstName = x$pt.method,
CIEstName = x$ci.method,
Alpha = x$alpha) # c() just adds to the list x
out <- lapply(attr(object,"fullList"), sumf)
#attributes(out) <- list(class = "summary.pooledBinList", names = attr(object,"names"),group.var = attr(object,"group.var"))
#attributes(out) <- list(class = "summary.pooledBinList",
# names = grp.names,group.var = attr(object,"group.var"))
n <- length(out)
group.dat <- as.data.frame(attr(object,"group.dat"))
if(ncol(group.dat)==0){
group.dat <- data.frame(Group=0)
} else {
names(group.dat) <- names(object)[1:length(attr(object,"group.var"))]
}
ans <- cbind(group.dat,
data.frame(
P = rep(0,n),
Lower = rep(0,n),
Upper = rep(0,n),
Scale = rep(1,n),
N = rep(0,n),
NumPools = rep(0,n),
NumPosPools = rep(0,n))
)
# ans <- data.frame(Group = grp.names,
# P = rep(0,n),
# Lower = rep(0,n),
# Upper = rep(0,n),
# Scale = rep(1,n),
# N = rep(0,n),
# NumPools = rep(0,n),
# NumPosPools = rep(0,n))
#if(!is.null(attr(object,"group.var"))) names(ans)[1] <- attr(object,"group.var")
#if(attr(object,"group.var") != "") names(ans)[1:length(attr(object,"group.var"))] <- attr(object,"group.var")
for(i in 1:n)
ans[i,(length(attr(object,"group.var"))-1) + (2:8)] <- c(out[[i]]$scale * c(out[[i]]$p, out[[i]]$lcl,out[[i]]$ucl), out[[i]]$scale,
out[[i]]$N, out[[i]]$NumPools, out[[i]]$NumPosPools)
if(attr(object,"n.groups") == 1) ans$Group <- NULL
if(simple) return(ans)
structure(ans, class = "summary.pIR",
df = as.data.frame(unclass(ans)),
fullList = attr(object,"fullList"),
#names = grp.names,
group.var = attr(object,"group.var"),
call = attr(object,"call"))
}
"print.summary.pIR" <- function(x, ...){
cat("\nEstimation of Binomial Proportion for Pooled Data\n")
cat(paste0("\nCall: ", deparse(attr(x,"call"),width.cutoff = 100),"\n\n"))
switch(attr(x,"fullList")[[1]]$pt.method,
"firth" = PtEstName <- "Firth's Correction",
"gart" = PtEstName <- "Gart's Correction",
"bc-mle" = PtEstName <- "Gart's Correction",
"mle" = PtEstName <- "Maximum Likelihood",
"mir" = PtEstName <- "Minimum Infection Rate"
)
switch(attr(x,"fullList")[[1]]$ci.method,
"skew-score" = CIEstName <- "Skew-Corrected Score (Gart)",
"bc-skew-score" = CIEstName <- "Bias- & Skew-Corrected Score (Gart)",
"score" = CIEstName <- "Score",
"lrt" = CIEstName <- "Likelihood Ratio Test Inversion",
"wald" = CIEstName <- "Wald",
"mir" = CIEstName <- "Minimum Infection Rate"
)
cat(paste("Point estimator :",PtEstName,"\n"))
cat(paste("CI method :",CIEstName,"\n"))
cat(paste("Confidence coefficient : ",100*(1-attr(x,"fullList")[[1]]$alpha),"%\n\n",sep=""))
print(attr(x,"df"))
invisible(x)
}
"plot.pIR" <-
function(x,pch=16,refline=TRUE,printR2=TRUE,layout=NULL,...){
x.lst <- attr(x,"fullList")
n.groups <- length(x.lst)
n.groups.root <- ceiling(sqrt(n.groups))
if(is.null(layout)) layout <- c(n.groups.root, n.groups.root)
par(mfrow = layout)
for(i in 1:n.groups){
# reference: Chen & Swallow
if(all(x.lst[[i]]$n==1)) {
xmn <- as.list(by(x.lst[[i]]$x, x.lst[[i]]$m, function(x) c(sum(x),length(x))))
m <- as.numeric(names(xmn))
xx <- sapply(xmn,function(x) x[1])
n <- sapply(xmn,function(x) x[2])
} else {
xx <- x.lst[[i]]$x
m <- x.lst[[i]]$m
n <- x.lst[[i]]$n
}
y <- log((xx+0.5)/(n+0.5))
cc <- lm(y ~ m)
plot(m,y,xlab="Pool Size",ylab="log((x+0.5)/(n+0.5))",...)
if(refline) abline(cc)
#if(printR2) title(bquote(.(names(x)[[i]]) ~ ":" ~ R^2 ~ "=" ~ .(round(summary(cc)$r.squared,2))), cex = 0.75)
#if(printR2) title(bquote(.(attr(x,"group.names")[[i]]) ~ ":" ~ R^2 ~ "=" ~ .(round(summary(cc)$r.squared,2))), cex = 0.75)
#else title(names(x)[[i]],cex=0.75)
if(n.groups > 1){
title.group <- as.character(attr(x,"group.names"))[[i]]
if(printR2) title(bquote(.(title.group) ~ ":" ~ R^2 ~ "=" ~ .(round(summary(cc)$r.squared,2))), cex = 0.75)
else title(bquote(.(title.group)), cex = 0.75)
} else {
if(printR2) title(bquote(R^2 ~ "=" ~ .(round(summary(cc)$r.squared,2))), cex = 0.75)
else title("Model Diagnostic",cex=0.75)
}
}
invisible(x)
}
"as.data.frame.pIR" <- function(x){
as.data.frame(unclass(x))
}
bjbiggerstaff/PooledInfRate documentation built on Jan. 19, 2024, 6:54 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
"pIR" <- function(x,...){
UseMethod("pIR")
}
"pIR.default" <- function(x,m,n=rep(1,length(x)), group,
pt.method = c("firth","gart","bc-mle","mle","mir"),
ci.method = c("skew-score","bc-skew-score","score","lrt","wald","mir"),
scale=1, alpha=0.05, tol=.Machine$double.eps^0.5, ...) {
call <- match.call()
call[[1]] <- as.name("pIR")
xmn.nomiss <- (!is.na(x) & !is.na(m) & !is.na(n))
x <- x[xmn.nomiss]
m <- m[xmn.nomiss]
n <- n[xmn.nomiss]
mn.nozero <- (m>0 & n>0)
x <- x[mn.nozero]
m <- m[mn.nozero]
n <- n[mn.nozero]
if(!missing(group)){
group.var <- "Group" #deparse(substitute(group))
xmng.nomiss <- (!is.na(x) & !is.na(m) & !is.na(n) & !is.na(group))
group <- group[xmng.nomiss]
groups <- unique(group)
nGroups <- length(groups)
group.dat <- data.frame(Group=group)
} else {
xmng.nomiss <- (!is.na(x) & !is.na(m) & !is.na(n))
group <- rep("SINGLE",length(x[xmng.nomiss]))
groups <- unique(group)
nGroups <- length(groups)
group.var <- ""
group.dat <- data.frame(Group=rep(1,length(x[xmng.nomiss])))
}
ans <- vector(mode="list",length=nGroups)
for(i in 1:nGroups){
ans[[i]] <- pooledBin.fit(x[group==groups[i]],
m[group==groups[i]],
n[group==groups[i]],
pt.method=pt.method,
ci.method=ci.method,
scale=scale,alpha=alpha,tol=tol)
}
ans.lst <- ans
ans <- cbind(group.dat[!duplicated(group),],
data.frame(
P = rep(0, nGroups),
Lower = rep(0, nGroups),
Upper = rep(0, nGroups),
Scale = rep(1, nGroups))
)
#ans <- data.frame(Group = groups,
# P = rep(0, nGroups),
# Lower = rep(0, nGroups),
# Upper = rep(0, nGroups),
# Scale = rep(1, nGroups))
names(ans)[1] <- ifelse(group.var != "",group.var,"Group")
#if (group.var != "") names(ans)[1] <- group.var
for (i in 1:nGroups) ans[i, 2:5] <- ans.lst[[i]]$scale * c(ans.lst[[i]]$p, ans.lst[[i]]$lcl, ans.lst[[i]]$ucl, 1)
if (all(ans$Scale == 1)) ans$Scale <- NULL
if(nGroups == 1) ans$Group <- NULL
# fix row names after subsetting
group.dat <- as.data.frame(group.dat[!duplicated(group),])
rownames(group.dat) <- 1:nrow(group.dat)
ans <- structure(ans, class = "pooledBin",fullList = ans.lst,
group.names = groups, group.var = group.var,n.groups = nGroups,
group.dat = group.dat,
scale=scale,call=call)
ans
}
# "pIR.formula" <- function(x, data,
# pt.method = c("firth","gart","bc-mle","mle","mir"),
# ci.method = c("skew-score","bc-skew-score","score","lrt","wald","mir"),
# scale=1, alpha=0.05, tol=.Machine$double.eps^0.5, ...){
# call <- match.call()
# call[[1]] <- as.name("pIR")
#
# if(missing(data))
# data <- environment(x)
#
# vars <- pooledBinParseFormula(x, data)
# if(any(sapply(vars[1:3],length)>1)) stop("only variable names permitted in formula; perhaps use the default call")
# # omit records with missing data -- note, if data contains records missing
# # anywhere (even not in X, M, N, Group, they are omitted), so care should be
# # used in subsetting before the call to be assured of the desired analysis
# # restrict to only those variables needed before subsetting to avoid that issue.
# if(!missing(data)){
# # restrict the data to the variables needed before removing records with any NA
# # --this doesn't help when no data are specified, so that's dealt with below
# data <- data[,unlist(vars)]
# data <- na.omit(data)
# }
#
# # retrieve values from data using the character name
# # use the eval(parse(text=XX), data) construct in case data is left unspecified,
# # and because we have the names of the variables available
# # from pooledBinParseFormula
# x <- eval(parse(text=vars$x), data)
# m <- eval(parse(text=vars$m), data)
# if(!is.null(vars$n))
# n <- eval(parse(text=vars$n), data)
# else n <- rep(1,length(x)) # default n
#
# if(!is.null(vars$group[1])){
# # NEW as of 5/24/2022 to allow multiple grouping variables
# n.group.var <- length(vars$group)
# group.dat <- as.data.frame(matrix(nrow=length(x),ncol=n.group.var))
# names(group.dat) <- vars$group
#
# for(i in 1:n.group.var) group.dat[,i] <- eval(parse(text=vars$group[i]),data)
#
# #print(interaction(group.dat,drop=TRUE)) # drop unused levels
# #return(group.dat)
#
# # end of NEW
#
# #group <- eval(parse(text=vars$group), data)
# group <- interaction(group.dat,drop=TRUE)
# groups <- unique(group)
# nGroups <- length(groups)
# group.var <- vars$group
# } else {
# group <- rep("SINGLE",length(x))
# groups <- unique(group)
# nGroups <- 1
# group.var <- ""
# group.dat <- data.frame()
# }
#
# # restrict to data with no missing x, m, n, group
# if(missing(data)){
# mn.nozero <- (m>0 & n>0)
# x <- x[mn.nozero]
# m <- m[mn.nozero]
# n <- n[mn.nozero]
#
# if(nGroups == 1){
# xmng.nomiss <- (!is.na(x) & !is.na(m) & !is.na(n))
# } else {
# xmng.nomiss <- (!is.na(x) & !is.na(m) & !is.na(n) & !is.na(group))
# }
# x <- x[xmng.nomiss]
# m <- m[xmng.nomiss]
# n <- n[xmng.nomiss]
#
# group <- group[xmng.nomiss]
# groups <- unique(group)
# nGroups <- length(groups)
# }
#
# mn.nozero <- (m>0 & n>0)
# x <- x[mn.nozero]
# m <- m[mn.nozero]
# n <- n[mn.nozero]
# group <- group[mn.nozero]
# groups <- unique(group)
# nGroups <- length(groups)
#
# #if(nGroups > 1){
# ans <- vector(mode="list",length=nGroups)
# for(i in 1:nGroups){
# ans[[i]] <- pooledBin.fit(x[group==groups[i]],
# m[group==groups[i]],
# n[group==groups[i]],
# pt.method=pt.method,
# ci.method=ci.method,
# scale=scale,alpha=alpha,tol=tol)
# }
# names(ans) <- groups
# ans.lst <- ans
#
# # ans <- data.frame(Group = groups,
# # P = rep(0, nGroups),
# # Lower = rep(0, nGroups),
# # Upper = rep(0, nGroups),
# # Scale = rep(1, nGroups))
# ans <- cbind(group.dat[!duplicated(group),],
# data.frame(
# P = rep(0, nGroups),
# Lower = rep(0, nGroups),
# Upper = rep(0, nGroups),
# Scale = rep(1, nGroups))
# )
# if (!is.null(vars$group)){
# #names(ans)[1] <- paste0(vars$group,collapse=".") # new as of 5/26/2022 for multiple grouping variables
# names(ans)[1:length(vars$group)] <- vars$group
# }
# for (i in 1:nGroups) ans[i,(ncol(group.dat)-1)+ (2:5)] <- ans.lst[[i]]$scale * c(ans.lst[[i]]$p, ans.lst[[i]]$lcl, ans.lst[[i]]$ucl, 1)
# if (all(ans$Scale == 1)) ans$Scale <- NULL
#
# if(nGroups == 1) ans$Group <- NULL
#
# # fix row names after subsetting
# group.dat <- as.data.frame(group.dat[!duplicated(group),])
# rownames(group.dat) <- 1:nrow(group.dat)
#
# ans <- structure(ans, class = "pIR", fullList = ans.lst,
# x.var = vars$x, m.var = vars$m, n.var = vars$n,
# group.names = groups, group.var = group.var,
# group.dat = group.dat,
# n.groups = nGroups, scale=scale,call=call)
# ans
# }
"pIR.formula" <- function(x, data,
pt.method = c("firth","gart","bc-mle","mle","mir"),
ci.method = c("skew-score","bc-skew-score","score","lrt","wald","mir"),
scale=1, alpha=0.05, tol=.Machine$double.eps^0.5, ...){
call <- match.call()
call[[1]] <- as.name("pooledBin")
have.df <- (!missing(data))
if(!have.df){
data <- environment(as.formula(x))
}
#data <- environment(x)
vars <- pooledBinParseFormula(x, data)
if(any(sapply(vars[1:3],length)>1)) stop("only variable names permitted in formula; perhaps use the default call")
# omit records with missing data -- note, if data contains records missing
# anywhere (even not in X, M, N, Group, they are omitted), so care should be
# used in subsetting before the call to be assured of the desired analysis
# restrict to only those variables needed before subsetting to avoid that issue.
#if(!missing(data)){
if(have.df){
# restrict the data to the variables needed before removing records with any NA
# --this doesn't help when no data are specified, so that's dealt with below
data <- data[,unlist(vars)]
data <- na.omit(data)
}
# retrieve values from data using the character name
# use the eval(parse(text=XX), data) construct in case data is left unspecified,
# and because we have the names of the variables available
# from pooledBinParseFormula
x <- eval(parse(text=vars$x), data)
m <- eval(parse(text=vars$m), data)
if(!is.null(vars$n))
n <- eval(parse(text=vars$n), data)
else n <- rep(1,length(x)) # default n
if(!is.null(vars$group[1])){
# NEW as of 5/24/2022 to allow multiple grouping variables
n.group.var <- length(vars$group)
group.dat <- as.data.frame(matrix(nrow=length(x),ncol=n.group.var))
names(group.dat) <- vars$group
for(i in 1:n.group.var) group.dat[,i] <- eval(parse(text=vars$group[i]),data)
#print(interaction(group.dat,drop=TRUE)) # drop unused levels
#return(group.dat)
# end of NEW
#group <- eval(parse(text=vars$group), data)
group <- interaction(group.dat,drop=TRUE)
groups <- unique(group)
nGroups <- length(groups)
group.var <- vars$group
} else {
group <- rep("SINGLE",length(x))
groups <- unique(group)
nGroups <- 1
group.var <- ""
group.dat <- data.frame()
}
# restrict to data with no missing x, m, n, group
#if(missing(data)){
mn.nozero <- (m>0 & n>0)
x <- x[mn.nozero]
m <- m[mn.nozero]
n <- n[mn.nozero]
if(nGroups == 1){
xmng.nomiss <- (!is.na(x) & !is.na(m) & !is.na(n))
} else {
xmng.nomiss <- (!is.na(x) & !is.na(m) & !is.na(n) & !is.na(group))
}
x <- x[xmng.nomiss]
m <- m[xmng.nomiss]
n <- n[xmng.nomiss]
group <- group[xmng.nomiss]
groups <- unique(group)
nGroups <- length(groups)
#}
# mn.nozero <- (m>0 & n>0)
# x <- x[mn.nozero]
# m <- m[mn.nozero]
# n <- n[mn.nozero]
# group <- group[mn.nozero]
# groups <- unique(group)
# nGroups <- length(groups)
#if(nGroups > 1){
ans <- vector(mode="list",length=nGroups)
for(i in 1:nGroups){
ans[[i]] <- pooledBin.fit(x[group==groups[i]],
m[group==groups[i]],
n[group==groups[i]],
pt.method=pt.method,
ci.method=ci.method,
scale=scale,alpha=alpha,tol=tol)
}
names(ans) <- groups
ans.lst <- ans
# ans <- data.frame(Group = groups,
# P = rep(0, nGroups),
# Lower = rep(0, nGroups),
# Upper = rep(0, nGroups),
# Scale = rep(1, nGroups))
ans <- cbind(group.dat[!duplicated(group),],
data.frame(
P = rep(0, nGroups),
Lower = rep(0, nGroups),
Upper = rep(0, nGroups),
Scale = rep(1, nGroups))
)
if (!is.null(vars$group)){
#names(ans)[1] <- paste0(vars$group,collapse=".") # new as of 5/26/2022 for multiple grouping variables
names(ans)[1:length(vars$group)] <- vars$group
}
for (i in 1:nGroups) ans[i,(ncol(group.dat)-1)+ (2:5)] <- ans.lst[[i]]$scale * c(ans.lst[[i]]$p, ans.lst[[i]]$lcl, ans.lst[[i]]$ucl, 1)
if (all(ans$Scale == 1)) ans$Scale <- NULL
if(nGroups == 1) ans$Group <- NULL
# fix row names after subsetting
group.dat <- as.data.frame(group.dat[!duplicated(group),])
rownames(group.dat) <- 1:nrow(group.dat)
ans <- structure(ans, class = "pIR", fullList = ans.lst,
x.var = vars$x, m.var = vars$m, n.var = vars$n,
group.names = groups, group.var = group.var,
group.dat = group.dat,
n.groups = nGroups, scale=scale,call=call)
ans
}
"print.pIR" <- function(x, ...){
print(as.data.frame(unclass(x)),...)
invisible(x)
}
# to build the package, had to use a version of "[.data.frame" without the .Internal(copyDFattr(xx,x)) call
# this just copied x as a list along with the attributes of the data frame
#"[.pIR" <- subsetDF #get("[.data.frame")
#"[.summary.pIR" <- subsetDF #get("[.data.frame")
"[.pIR" <- function (x, i, j, drop = if (missing(i)) TRUE else length(cols) == 1) {
mdrop <- missing(drop)
Narg <- nargs() - !mdrop
has.j <- !missing(j)
if (!all(names(sys.call()) %in% c("", "drop")) &&
!isS4(x))
warning("named arguments other than 'drop' are discouraged")
if (Narg < 3L) {
if (!mdrop)
warning("'drop' argument will be ignored")
if (missing(i))
return(x)
if (is.matrix(i))
return(as.matrix(x)[i])
nm <- names(x)
if (is.null(nm))
nm <- character()
if (!is.character(i) && anyNA(nm)) {
names(nm) <- names(x) <- seq_along(x)
y <- NextMethod("[")
cols <- names(y)
if (anyNA(cols))
stop("undefined columns selected")
cols <- names(y) <- nm[cols]
}
else {
y <- NextMethod("[")
cols <- names(y)
if (!is.null(cols) && anyNA(cols))
stop("undefined columns selected")
}
if (anyDuplicated(cols))
names(y) <- make.unique(cols)
attr(y, "row.names") <- .row_names_info(x, 0L)
attr(y, "class") <- oldClass(x)
return(y)
}
if (missing(i)) {
if (drop && !has.j && length(x) == 1L)
return(.subset2(x, 1L))
nm <- names(x)
if (is.null(nm))
nm <- character()
if (has.j && !is.character(j) && anyNA(nm)) {
names(nm) <- names(x) <- seq_along(x)
y <- .subset(x, j)
cols <- names(y)
if (anyNA(cols))
stop("undefined columns selected")
cols <- names(y) <- nm[cols]
}
else {
y <- if (has.j)
.subset(x, j)
else x
cols <- names(y)
if (anyNA(cols))
stop("undefined columns selected")
}
if (drop && length(y) == 1L)
return(.subset2(y, 1L))
if (anyDuplicated(cols))
names(y) <- make.unique(cols)
nrow <- .row_names_info(x, 2L)
if (drop && !mdrop && nrow == 1L)
return(structure(y, class = NULL, row.names = NULL))
else {
attr(y, "class") <- oldClass(x)
attr(y, "row.names") <- .row_names_info(x,
0L)
return(y)
}
}
xx <- x
cols <- names(xx)
x <- vector("list", length(x))
#x <- .Internal(copyDFattr(xx, x))
x <- as.list(xx)
attributes(x) <- attributes(xx)
oldClass(x) <- attr(x, "row.names") <- NULL
if (has.j) {
nm <- names(x)
if (is.null(nm))
nm <- character()
if (!is.character(j) && anyNA(nm))
names(nm) <- names(x) <- seq_along(x)
x <- x[j]
cols <- names(x)
if (drop && length(x) == 1L) {
if (is.character(i)) {
rows <- attr(xx, "row.names")
i <- pmatch(i, rows, duplicates.ok = TRUE)
}
xj <- .subset2(.subset(xx, j), 1L)
return(if (length(dim(xj)) != 2L) xj[i] else xj[i,
, drop = FALSE])
}
if (anyNA(cols))
stop("undefined columns selected")
if (!is.null(names(nm)))
cols <- names(x) <- nm[cols]
nxx <- structure(seq_along(xx), names = names(xx))
sxx <- match(nxx[j], seq_along(xx))
}
else sxx <- seq_along(x)
rows <- NULL
if (is.character(i)) {
rows <- attr(xx, "row.names")
i <- pmatch(i, rows, duplicates.ok = TRUE)
}
for (j in seq_along(x)) {
xj <- xx[[sxx[j]]]
x[[j]] <- if (length(dim(xj)) != 2L)
xj[i]
else xj[i, , drop = FALSE]
}
if (drop) {
n <- length(x)
if (n == 1L)
return(x[[1L]])
if (n > 1L) {
xj <- x[[1L]]
nrow <- if (length(dim(xj)) == 2L)
dim(xj)[1L]
else length(xj)
drop <- !mdrop && nrow == 1L
}
else drop <- FALSE
}
if (!drop) {
if (is.null(rows))
rows <- attr(xx, "row.names")
rows <- rows[i]
if ((ina <- anyNA(rows)) | (dup <- anyDuplicated(rows))) {
if (!dup && is.character(rows))
dup <- "NA" %in% rows
if (ina)
rows[is.na(rows)] <- "NA"
if (dup)
rows <- make.unique(as.character(rows))
}
if (has.j && anyDuplicated(nm <- names(x)))
names(x) <- make.unique(nm)
if (is.null(rows))
rows <- attr(xx, "row.names")[i]
attr(x, "row.names") <- rows
oldClass(x) <- oldClass(xx)
}
x
}
"[.summary.pIR" <- function (x, i, j, drop = if (missing(i)) TRUE else length(cols) == 1) {
mdrop <- missing(drop)
Narg <- nargs() - !mdrop
has.j <- !missing(j)
if (!all(names(sys.call()) %in% c("", "drop")) &&
!isS4(x))
warning("named arguments other than 'drop' are discouraged")
if (Narg < 3L) {
if (!mdrop)
warning("'drop' argument will be ignored")
if (missing(i))
return(x)
if (is.matrix(i))
return(as.matrix(x)[i])
nm <- names(x)
if (is.null(nm))
nm <- character()
if (!is.character(i) && anyNA(nm)) {
names(nm) <- names(x) <- seq_along(x)
y <- NextMethod("[")
cols <- names(y)
if (anyNA(cols))
stop("undefined columns selected")
cols <- names(y) <- nm[cols]
}
else {
y <- NextMethod("[")
cols <- names(y)
if (!is.null(cols) && anyNA(cols))
stop("undefined columns selected")
}
if (anyDuplicated(cols))
names(y) <- make.unique(cols)
attr(y, "row.names") <- .row_names_info(x, 0L)
attr(y, "class") <- oldClass(x)
return(y)
}
if (missing(i)) {
if (drop && !has.j && length(x) == 1L)
return(.subset2(x, 1L))
nm <- names(x)
if (is.null(nm))
nm <- character()
if (has.j && !is.character(j) && anyNA(nm)) {
names(nm) <- names(x) <- seq_along(x)
y <- .subset(x, j)
cols <- names(y)
if (anyNA(cols))
stop("undefined columns selected")
cols <- names(y) <- nm[cols]
}
else {
y <- if (has.j)
.subset(x, j)
else x
cols <- names(y)
if (anyNA(cols))
stop("undefined columns selected")
}
if (drop && length(y) == 1L)
return(.subset2(y, 1L))
if (anyDuplicated(cols))
names(y) <- make.unique(cols)
nrow <- .row_names_info(x, 2L)
if (drop && !mdrop && nrow == 1L)
return(structure(y, class = NULL, row.names = NULL))
else {
attr(y, "class") <- oldClass(x)
attr(y, "row.names") <- .row_names_info(x,
0L)
return(y)
}
}
xx <- x
cols <- names(xx)
x <- vector("list", length(x))
#x <- .Internal(copyDFattr(xx, x))
x <- as.list(xx)
attributes(x) <- attributes(xx)
oldClass(x) <- attr(x, "row.names") <- NULL
if (has.j) {
nm <- names(x)
if (is.null(nm))
nm <- character()
if (!is.character(j) && anyNA(nm))
names(nm) <- names(x) <- seq_along(x)
x <- x[j]
cols <- names(x)
if (drop && length(x) == 1L) {
if (is.character(i)) {
rows <- attr(xx, "row.names")
i <- pmatch(i, rows, duplicates.ok = TRUE)
}
xj <- .subset2(.subset(xx, j), 1L)
return(if (length(dim(xj)) != 2L) xj[i] else xj[i,
, drop = FALSE])
}
if (anyNA(cols))
stop("undefined columns selected")
if (!is.null(names(nm)))
cols <- names(x) <- nm[cols]
nxx <- structure(seq_along(xx), names = names(xx))
sxx <- match(nxx[j], seq_along(xx))
}
else sxx <- seq_along(x)
rows <- NULL
if (is.character(i)) {
rows <- attr(xx, "row.names")
i <- pmatch(i, rows, duplicates.ok = TRUE)
}
for (j in seq_along(x)) {
xj <- xx[[sxx[j]]]
x[[j]] <- if (length(dim(xj)) != 2L)
xj[i]
else xj[i, , drop = FALSE]
}
if (drop) {
n <- length(x)
if (n == 1L)
return(x[[1L]])
if (n > 1L) {
xj <- x[[1L]]
nrow <- if (length(dim(xj)) == 2L)
dim(xj)[1L]
else length(xj)
drop <- !mdrop && nrow == 1L
}
else drop <- FALSE
}
if (!drop) {
if (is.null(rows))
rows <- attr(xx, "row.names")
rows <- rows[i]
if ((ina <- anyNA(rows)) | (dup <- anyDuplicated(rows))) {
if (!dup && is.character(rows))
dup <- "NA" %in% rows
if (ina)
rows[is.na(rows)] <- "NA"
if (dup)
rows <- make.unique(as.character(rows))
}
if (has.j && anyDuplicated(nm <- names(x)))
names(x) <- make.unique(nm)
if (is.null(rows))
rows <- attr(xx, "row.names")[i]
attr(x, "row.names") <- rows
oldClass(x) <- oldClass(xx)
}
x
}
"summary.pIR" <- function(object, simple = FALSE, ...){
grp.names <- as.character(attr(object,"group.names"))
"sumf" <- function(x) c(x, N = sum(x$n * x$m), NumPools = sum(x$n), NumPosPools = sum(x$x),
PtEstName = x$pt.method,
CIEstName = x$ci.method,
Alpha = x$alpha) # c() just adds to the list x
out <- lapply(attr(object,"fullList"), sumf)
#attributes(out) <- list(class = "summary.pooledBinList", names = attr(object,"names"),group.var = attr(object,"group.var"))
#attributes(out) <- list(class = "summary.pooledBinList",
# names = grp.names,group.var = attr(object,"group.var"))
n <- length(out)
group.dat <- as.data.frame(attr(object,"group.dat"))
if(ncol(group.dat)==0){
group.dat <- data.frame(Group=0)
} else {
names(group.dat) <- names(object)[1:length(attr(object,"group.var"))]
}
ans <- cbind(group.dat,
data.frame(
P = rep(0,n),
Lower = rep(0,n),
Upper = rep(0,n),
Scale = rep(1,n),
N = rep(0,n),
NumPools = rep(0,n),
NumPosPools = rep(0,n))
)
# ans <- data.frame(Group = grp.names,
# P = rep(0,n),
# Lower = rep(0,n),
# Upper = rep(0,n),
# Scale = rep(1,n),
# N = rep(0,n),
# NumPools = rep(0,n),
# NumPosPools = rep(0,n))
#if(!is.null(attr(object,"group.var"))) names(ans)[1] <- attr(object,"group.var")
#if(attr(object,"group.var") != "") names(ans)[1:length(attr(object,"group.var"))] <- attr(object,"group.var")
for(i in 1:n)
ans[i,(length(attr(object,"group.var"))-1) + (2:8)] <- c(out[[i]]$scale * c(out[[i]]$p, out[[i]]$lcl,out[[i]]$ucl), out[[i]]$scale,
out[[i]]$N, out[[i]]$NumPools, out[[i]]$NumPosPools)
if(attr(object,"n.groups") == 1) ans$Group <- NULL
if(simple) return(ans)
structure(ans, class = "summary.pIR",
df = as.data.frame(unclass(ans)),
fullList = attr(object,"fullList"),
#names = grp.names,
group.var = attr(object,"group.var"),
call = attr(object,"call"))
}
"print.summary.pIR" <- function(x, ...){
cat("\nEstimation of Binomial Proportion for Pooled Data\n")
cat(paste0("\nCall: ", deparse(attr(x,"call"),width.cutoff = 100),"\n\n"))
switch(attr(x,"fullList")[[1]]$pt.method,
"firth" = PtEstName <- "Firth's Correction",
"gart" = PtEstName <- "Gart's Correction",
"bc-mle" = PtEstName <- "Gart's Correction",
"mle" = PtEstName <- "Maximum Likelihood",
"mir" = PtEstName <- "Minimum Infection Rate"
)
switch(attr(x,"fullList")[[1]]$ci.method,
"skew-score" = CIEstName <- "Skew-Corrected Score (Gart)",
"bc-skew-score" = CIEstName <- "Bias- & Skew-Corrected Score (Gart)",
"score" = CIEstName <- "Score",
"lrt" = CIEstName <- "Likelihood Ratio Test Inversion",
"wald" = CIEstName <- "Wald",
"mir" = CIEstName <- "Minimum Infection Rate"
)
cat(paste("Point estimator :",PtEstName,"\n"))
cat(paste("CI method :",CIEstName,"\n"))
cat(paste("Confidence coefficient : ",100*(1-attr(x,"fullList")[[1]]$alpha),"%\n\n",sep=""))
print(attr(x,"df"))
invisible(x)
}
"plot.pIR" <-
function(x,pch=16,refline=TRUE,printR2=TRUE,layout=NULL,...){
x.lst <- attr(x,"fullList")
n.groups <- length(x.lst)
n.groups.root <- ceiling(sqrt(n.groups))
if(is.null(layout)) layout <- c(n.groups.root, n.groups.root)
par(mfrow = layout)
for(i in 1:n.groups){
# reference: Chen & Swallow
if(all(x.lst[[i]]$n==1)) {
xmn <- as.list(by(x.lst[[i]]$x, x.lst[[i]]$m, function(x) c(sum(x),length(x))))
m <- as.numeric(names(xmn))
xx <- sapply(xmn,function(x) x[1])
n <- sapply(xmn,function(x) x[2])
} else {
xx <- x.lst[[i]]$x
m <- x.lst[[i]]$m
n <- x.lst[[i]]$n
}
y <- log((xx+0.5)/(n+0.5))
cc <- lm(y ~ m)
plot(m,y,xlab="Pool Size",ylab="log((x+0.5)/(n+0.5))",...)
if(refline) abline(cc)
#if(printR2) title(bquote(.(names(x)[[i]]) ~ ":" ~ R^2 ~ "=" ~ .(round(summary(cc)$r.squared,2))), cex = 0.75)
#if(printR2) title(bquote(.(attr(x,"group.names")[[i]]) ~ ":" ~ R^2 ~ "=" ~ .(round(summary(cc)$r.squared,2))), cex = 0.75)
#else title(names(x)[[i]],cex=0.75)
if(n.groups > 1){
title.group <- as.character(attr(x,"group.names"))[[i]]
if(printR2) title(bquote(.(title.group) ~ ":" ~ R^2 ~ "=" ~ .(round(summary(cc)$r.squared,2))), cex = 0.75)
else title(bquote(.(title.group)), cex = 0.75)
} else {
if(printR2) title(bquote(R^2 ~ "=" ~ .(round(summary(cc)$r.squared,2))), cex = 0.75)
else title("Model Diagnostic",cex=0.75)
}
}
invisible(x)
}
"as.data.frame.pIR" <- function(x){
as.data.frame(unclass(x))
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.