binomInference.exactLR <- function(y, n, null.hypothesis=NA, test.type="two.sided", conf.level=0.95, version=F) { vrsn <- "20121107" if (version) return(vrsn) alpha <- (1 - conf.level) / 2 if(y==0) elo <-0 else { elo <- qbeta(alpha, y, n-y+1) } if(y==n) ehi <- 1 else { ehi <- qbeta(1-alpha,y+1,n-y) } if(!is.na(null.hypothesis)) { plo <- pbinom(y, n, null.hypothesis) phi <- 1 - plo + dbinom(y, n, null.hypothesis) pdf <- dbinom(0:n, n, null.hypothesis) if (plo < phi) { u <- (0:n) >= (n * null.hypothesis) p2 <- plo } else { u <- (0:n) <= (n * null.hypothesis) p2 <- phi } u <- u & (pdf <= dbinom(y, n, null.hypothesis) * (1 + 1e-7)) if(sum(u) > 0) p2 <- min(p2 + sum(pdf[u]),1) if(test.type=="less") p.value <- plo else if(test.type=="greater") p.value <- phi else p.value <- p2 } else p.value <- NA rslt <- c(n, y/n, NA, NA, elo, ehi, p.value) names(rslt) <- c("n", "Est", "SE", "Statistic", "CIlo", "CIhi", "p.value") attr(rslt,"method") = "Exact(LR)" rslt } binomInference.exactTail <- function(y, n, null.hypothesis=NA, test.type="two.sided", conf.level=0.95, version=F) { vrsn <- "20121107" if (version) return(vrsn) alpha <- (1 - conf.level) / 2 if(y==0) elo <-0 else { elo <- qbeta(alpha, y, n-y+1) } if(y==n) ehi <- 1 else { ehi <- qbeta(1-alpha,y+1,n-y) } if(!is.na(null.hypothesis)) { plo <- pbinom(y, n, null.hypothesis) phi <- 1 - plo + dbinom(y, n, null.hypothesis) cdf <- pbinom(0:n, n, null.hypothesis) if (plo < phi) { cdf <- 1 - cdf p2 <- plo } else p2 <- phi u <- cdf <= p2 * (1 + 1e-7) if(sum(u) > 0) p2 <- min(p2 + max(cdf[u]),1) if(test.type=="less") p.value <- plo else if(test.type=="greater") p.value <- phi else p.value <- p2 } else p.value <- NA rslt <- c(n, y/n, NA, NA, elo, ehi, p.value) names(rslt) <- c("n", "Est", "SE", "Statistic", "CIlo", "CIhi", "p.value") attr(rslt,"method") = "Exact(tail)" rslt } binomInference.halfP <- function(y, n, null.hypothesis=NA, test.type="two.sided", conf.level=0.95, version=F) { vrsn <- "20121107" if (version) return(vrsn) alpha <- (1 - conf.level) / 2 if(y==0) hlo <- 0 else { f <- function(x, n, y, alpha) pbinom(n-y, n, 1-x) - dbinom(n-y, n, 1-x) / 2 - alpha hlo <- uniroot(f, c(0,1), n=n, y=y, alpha=alpha, tol=1e-10)$root } if(y==n) hhi <- 1 else { f <- function (x, n, y, alpha) pbinom(y, n, x) - dbinom(y, n, x) / 2- alpha hhi <- uniroot (f, c(0,1), n=n, y=y, alpha=alpha, tol=1e-10)$root } if(!is.na(null.hypothesis)) { plo <- pbinom(y, n, null.hypothesis) - dbinom(y, n, null.hypothesis)/2 phi <- 1 - plo + dbinom(y, n, null.hypothesis) / 2 cdf <- pbinom(0:n, n, null.hypothesis) - dbinom(0:n, n, null.hypothesis) / 2 if (plo < phi) { cdf <- 1 - cdf p2 <- plo } else p2 <- phi u <- cdf <= p2 * (1 + 1e-7) if(sum(u) > 0) p2 <- min(p2 + max(cdf[u]),1) if(test.type=="less") p.value <- plo else if(test.type=="greater") p.value <- phi else p.value <- p2 } else p.value <- NA rslt <- c(n, y/n, NA, NA, elo, ehi, p.value) names(rslt) <- c("n", "Est", "SE", "Statistic", "CIlo", "CIhi", "p.value") attr(rslt,"method") = "HalfP" rslt } binomInference.jeffreys <- function(y, n, null.hypothesis=NA, test.type="two.sided", conf.level=0.95, version=F) { vrsn <- "20121107" if (version) return(vrsn) alpha <- (1 - conf.level) / 2 phat <- qbeta(0.5, y+0.5, n-y+0.5) if(y==0) elo <-0 else { elo <- qbeta(alpha, y+0.5, n-y+0.5) } if(y==n) ehi <- 1 else { ehi <- qbeta(1-alpha,y+0.5,n-y+0.5) } if(!is.na(null.hypothesis)) { plo <- pbeta(null.hypothesis, y+0.5, n-y+0.5) phi <- 1 - plo p2 <- 2 * min(plo,phi) if(test.type=="less") p.value <- plo else if(test.type=="greater") p.value <- phi else p.value <- p2 } else p.value <- NA rslt <- c(n, phat, NA, NA, elo, ehi, p.value) names(rslt) <- c("n", "Est", "SE", "Statistic", "CIlo", "CIhi", "p.value") attr(rslt,"method") = "Jeffreys" rslt } binomInference.wald <- function(y, n, null.hypothesis=NA, test.type="two.sided", conf.level=0.95, version=F) { vrsn <- "20121107" if (version) return(vrsn) alpha <- (1 - conf.level) / 2 phat <- y / n wlo <- phat + qnorm(alpha)*sqrt(phat * (1-phat) / n) whi <- phat - qnorm(alpha)*sqrt(phat * (1-phat) / n) if(!is.na(null.hypothesis)) { SE <- sqrt(phat * (1 - phat) / n) Z <- (phat - null.hypothesis) / SE plo <- pnorm(Z) phi <- 1 - pnorm(Z) p2 <- 2 * min(plo,phi) if(test.type=="less") p.value <- plo else if(test.type=="greater") p.value <- phi else p.value <- p2 } else SE <- Z <- p.value <- NA rslt <- c(n, phat, SE, Z, wlo, whi, p.value) names(rslt) <- c("n", "Est", "SE", "Statistic", "CIlo", "CIhi", "p.value") attr(rslt,"method") = "Wald" rslt } binomInference.cwald <- function(y, n, null.hypothesis=NA, test.type="two.sided", conf.level=0.95, version=F) { vrsn <- "20121107" if (version) return(vrsn) alpha <- (1 - conf.level) / 2 phat <- y / n wlo <- phat - 1 / (2 * n) + qnorm(alpha)*sqrt(phat * (1-phat) / n) whi <- phat + 1 / (2 * n) - qnorm(alpha)*sqrt(phat * (1-phat) / n) if(!is.na(null.hypothesis)) { SE <- sqrt(phat * (1 - phat) / n) plo <- pnorm((phat + 1 / (2 * n) - null.hypothesis) / SE) phi <- 1 - pnorm((phat - 1 / (2 * n) - null.hypothesis) / SE) p2 <- 2 * min(plo,phi) if(test.type=="less") p.value <- plo else if(test.type=="greater") p.value <- phi else p.value <- p2 } else SE <- p.value <- NA rslt <- c(n, phat, SE, NA, wlo, whi, p.value) names(rslt) <- c("n", "Est", "SE", "Statistic", "CIlo", "CIhi", "p.value") attr(rslt,"method") = "cWald" rslt } binomInference.score <- function(y, n, null.hypothesis=NA, test.type="two.sided", conf.level=0.95, version=F) { vrsn <- "20121107" if (version) return(vrsn) alpha <- (1 - conf.level) / 2 phat <- y / n if(y==0) slo <- 0 else slo <- (2 * n * phat + qnorm(alpha)^2 - sqrt(4 * n * phat * (1 - phat) * qnorm(alpha)^2 + qnorm(alpha)^4)) / (2 * (n + qnorm(alpha)^2)) if(y==n) shi <- 1 else shi <- (2 * n * phat + qnorm(alpha)^2 + sqrt(4 * n * phat * (1 - phat) * qnorm(alpha)^2 + qnorm(alpha)^4)) / (2 * (n + qnorm(alpha)^2)) if(!is.na(null.hypothesis)) { SE <- sqrt(null.hypothesis * (1 - null.hypothesis) / n) Z <- (phat - null.hypothesis) / SE plo <- pnorm((phat - null.hypothesis) / SE) phi <- 1 - pnorm((phat - null.hypothesis) / SE) p2 <- 2 * min(plo,phi) if(test.type=="less") p.value <- plo else if(test.type=="greater") p.value <- phi else p.value <- p2 } else SE <- Z <- p.value <- NA rslt <- c(n, phat, SE, Z, slo, shi, p.value) names(rslt) <- c("n", "Est", "SE", "Statistic", "CIlo", "CIhi", "p.value") attr(rslt,"method") <- "Score" rslt } binomInference.cscore <- function(y, n, null.hypothesis=NA, test.type="two.sided", conf.level=0.95, version=F) { vrsn <- "20121107" if (version) return(vrsn) alpha <- (1 - conf.level) / 2 phat <- y / n phatL <- phat - 1/2/n phatU <- phat + 1/2/n if(y==0) slo <- 0 else slo <- (2 * n * phatL + qnorm(alpha)^2 - sqrt(4 * n * phatL * (1 - phatL) * qnorm(alpha)^2 + qnorm(alpha)^4)) / (2 * (n + qnorm(alpha)^2)) if(y==n) shi <- 1 else shi <- (2 * n * phatU + qnorm(alpha)^2 + sqrt(4 * n * phatU * (1 - phatU) * qnorm(alpha)^2 + qnorm(alpha)^4)) / (2 * (n + qnorm(alpha)^2)) if(!is.na(null.hypothesis)) { SE <- sqrt(null.hypothesis * (1 - null.hypothesis) / n) Z <- (phat - null.hypothesis) / SE plo <- pnorm((phat + 1 / (2 * n) - null.hypothesis) / SE) phi <- 1 - pnorm((phat - 1 / (2 * n) - null.hypothesis) / SE) p2 <- 2 * min(plo,phi) if(test.type=="less") p.value <- plo else if(test.type=="greater") p.value <- phi else p.value <- p2 } else SE <- p.value <- NA rslt <- c(n, phat, SE, NA, slo, shi, p.value) names(rslt) <- c("n", "Est", "SE", "Statistic", "CIlo", "CIhi", "p.value") rslt attr(rslt,"method") <- "cScore" rslt } binomInference.agresti <- function(y, n, null.hypothesis=NA, test.type="two.sided", conf.level=0.95, version=F) { vrsn <- "20121107" if (version) return(vrsn) alpha <- (1 - conf.level) / 2 y <- y + qnorm(alpha)^2 / 2 n <- n + qnorm(alpha)^2 phat <- y / n SE <- sqrt(phat * (1 - phat) / n) wlo <- phat + qnorm(alpha)*sqrt(phat * (1-phat) / n) whi <- phat - qnorm(alpha)*sqrt(phat * (1-phat) / n) if(!is.nA(null.hypothesis)) { Z <- (phat - null.hypothesis) / SE plo <- pnorm(Z) phi <- 1 - pnorm(Z) p2 <- 2 * min(plo,phi) if(test.type=="less") p.value <- plo else if(test.type=="greater") p.value <- phi else p.value <- p2 } else Z <- p.value <- NA rslt <- c(n, phat, SE, Z, wlo, whi, p.value) names(rslt) <- c("n", "Est", "SE", "Statistic", "CIlo", "CIhi", "p.value") attr(rslt,"method") = "Agresti" rslt } oneSample <- function (fnctl, # the distribution functional to be compared y, # the response variable in the regression analysis null.hypothesis=NA, test.type="two.sided", subset=rep(TRUE,N), conf.level=0.95, na.rm=TRUE, probs= 0.5, # for use with quantiles replaceZeroes=NULL, # for use with geometric means restriction=Inf, # for use with means, geometric means of censored data subjTime=rep(1,length(y)), # for use with rates method= NULL, above=NULL, below=NULL, labove=NULL, rbelow=NULL, interval=NULL, linterval=NULL, rinterval=NULL, lrinterval=NULL, # for use with proportions of discretized continuous data g1=1, g2=0, dispersion=1, # for use with method=="mean-variance" nbstrap=10000, resample="pairs", seed=0, # for use with method=="bootstrap" ..., version=FALSE) { # # Some notes on order of precedence: # If restriction is supplied, y is replaced by a restricted variable that is a "Surv" object only if is.Surv(y) # If any of the argurments above, ..., lrinterval are supplied and !is.Surv(y), y is replaced by an indicator variable # If more than one of the arguments above, ..., lrinterval are supplied, precedence is in the order given in the argument list (no warnings) # If none of the arguments above, ..., lrinterval are supplied and fnctl=="proportion" or fnctl=="odds", default is above=0 # Methods that are implemented depend on fnctl and type of data # If is.Surv(y) # If fnctl=="mean" or fnctl=="geometric mean" # "KM" (the default if is.Surv(y)) # The following methods can be used, but result in an error if any observations are censored prior to the restriction # "t.test" # "mean-variance" (which takes arguments g1, g2, and dispersion) # If fnctl=="median" or fnctl=="quantile" # "KM" (the default if is.Surv(y)) # "bootstrap" (which takes arguments nbstrap and resample=("pairs","independent")) # The following methods can be used, but result in an error if any observations are censored prior to the restriction # "sign" # If fnctl=="proportion" or fnctl=="odds" and # If fnctl=="mean" or fnctl=="geometric mean" # "t.test" (the default if !is.Surv(y)) # "mean-variance" (which takes arguments g1, g2, and dispersion) # "bootstrap" (which takes arguments nbstrap) # "KM" (the default if is.Surv(y)) # If fnctl=="median" or fnctl=="quantile" # "sign" (the default if !is.Surv(y)) # "bootstrap" (which takes arguments nbstrap) # "KM" (the default if is.Surv(y)) # If fnctl=="proportion" or fnctl=="odds" and # can be a character string or a list of character strings named $test and $interval; choices depend on fnctl # If fnctl=="mean" # method can be "t.test" (the default); "bootstrap" # # Output will include # - sample size # - functional name # - estimate # - confidence interval # - p value # - link function for analysis model # - estimate on analysis model scale # - standard error on analysis model scale (if relevant) vrsn <- "20121107" if (version) return(vrsn) chc.fnctl <- c("mean", "geometric mean", "proportion", "median", "quantile", "odds", "rate") name.fnctl <- c("Mean","GeomMn","Prop","Mdn",paste(format(100*probs[1]),"%ile",sep=""),"Odds","Rate") findx <- pmatch(fnctl, chc.fnctl) if (is.na(findx)) stop("unsupported functional") fnctl <- chc.fnctl[findx] isSurv <- is.Surv(y) isDate <- inherits(y,"Date") isEvent <- is.logical(y) || all(y[!is.na(y)] %in% 0:1) if(!isSurv && !isDate &&!isEvent && !is.numeric(y)) stop("y is an unsupported data type") if(isSurv) N <- dim(y)[1] else N <- length(y) if(restriction < Inf) { if(isSurv) { u <- y[,1] > restriction y[u,1] <- restriction y[u,2] <- 0 } else y[y > restriction] <- restriction rName <- paste("Summary measure restricted to", format(restriction)) } else rName <- NULL if (na.rm) u <- !is.na(y) else u <- rep(TRUE,N) ya <- y[subset & u] if(length(ya)==0) stop("no data to analyze") if(isSurv) n <- dim(ya)[1] else n <- length(ya) if(!na.rm) nMsng <- NA else nMsng <- sum(subset & !u) if(isSurv) default.method <- "KM" else default.method <- c("t.test", "t.test", "exact", "sign", "sign", "exact", "mean-variance")[findx] chc.method <- c("t.test", "exact", "exactLR", "exactTail", "wald", "cwald", "score", "cscore", "agresti", "jeffreys", "sign", "bootstrap", "mean-variance", "KM") if(is.null(method)) method <- default.method mindx <- pmatch(method, chc.method) if (is.na(mindx)) stop("unsupported method") method <- chc.method[mindx] thresholds <- NULL if (length(above)>0) thresholds <- rbind(thresholds,cbind(0,above,0,Inf)) if (length(below)>0) thresholds <- rbind(thresholds,cbind(0,-Inf,0,below)) if (length(labove)>0) thresholds <- rbind(thresholds,cbind(1,labove,0,Inf)) if (length(rbelow)>0) thresholds <- rbind(thresholds,cbind(0,-Inf,1,rbelow)) if (!is.null(interval)) { if (length(interval)==2) interval <- matrix(interval,ncol=2) if (dim(interval)[2]!=2) stop("intervals must be specified in a 2 column matrix") thresholds <- rbind(thresholds,cbind(0,interval[,1],0,interval[,2])) } if (!is.null(linterval)) { if (length(linterval)==2) linterval <- matrix(linterval,ncol=2) if (dim(linterval)[2]!=2) stop("intervals must be specified in a 2 column matrix") thresholds <- rbind(thresholds,cbind(1,linterval[,1],0,linterval[,2])) } if (!is.null(rinterval)) { if (length(rinterval)==2) rinterval <- matrix(rinterval,ncol=2) if (dim(rinterval)[2]!=2) stop("intervals must be specified in a 2 column matrix") thresholds <- rbind(thresholds,cbind(0,rinterval[,1],1,rinterval[,2])) } if (!is.null(lrinterval)) { if (length(lrinterval)==2) lrinterval <- matrix(lrinterval,ncol=2) if (dim(lrinterval)[2]!=2) stop("intervals must be specified in a 2 column matrix") thresholds <- rbind(thresholds,cbind(1,lrinterval[,1],1,lrinterval[,2])) } if (is.null(thresholds)) thresholds <- rbind(thresholds,cbind(0,0,0,Inf)) thresholds <- thresholds[1,,drop=F] if(fnctl %in% c("proportion","odds")) { if(isEvent) fName <- ifelse(fnctl=="proportion","Pr(Event)","Odds(Event)") else fName <- paste(sep="",ifelse(fnctl=="proportion","Pr","Odds"), ifelse(thresholds[,2] == -Inf, paste(sep="",ifelse(thresholds[,3]==0,"<","<="),format(thresholds[,4])), ifelse(thresholds[,4]==Inf, paste(sep="",ifelse(thresholds[,1]==0,">",">="),format(thresholds[,2])), paste(sep="", ifelse(thresholds[,1]==0,"(","["),format(thresholds[,2]),",",format(thresholds[,4]), ifelse(thresholds[,3]==0,")","]"))))) if (!isSurv) ya <- ifelse1(thresholds[,1]==0, ya > thresholds[,2], ya >= thresholds[,2]) & ifelse1(thresholds[,3]==0, ya < thresholds[,4], ya <= thresholds[,4]) } else fName <- name.fnctl[findx] chc.test <- c("greater","less","two.sided") tindx <- pmatch(test.type, chc.test) if (is.na(tindx)) stop("unsupported test type") test.type <- chc.test[tindx] if (is.na(null.hypothesis)) { hName <- NULL tindx <- 0 } else hName <- paste("Hypothesis test of", c("upper","lower","two-sided")[tindx], "alternative that", fName, c(">", "<", "<>")[tindx], format(null.hypothesis)) chc.resample <- c("pairs","independent") if(!is.null(resample)) { rindx <- pmatch(resample, chc.method) if (is.na(mindx)) stop("unsupported method") resample <- chc.resample[rindx] } nReplace <- NA if(isSurv) stop("KM based inference not yet implemented") else { if(fnctl=="mean") { if (method=="t.test") { mName <- "One sample t test" link <- "identity" etaName <- fName etaHat <- mean(ya) etaHatSE <- sqrt(var(ya) / n) etaNull <- null.hypothesis statistic <- (etaHat - etaNull) / etaHatSE df <- n - 1 if(test.type=="less") p.value <- pt(statistic,df) else if(test.type=="greater") p.value <- 1 - pt(statistic,df) else if(test.type=="two.sided") p.value <- 2 * pt(-abs(statistic),df) else p.value <- NA etaCIlo <- qt((1-conf.level)/2,df) * etaHatSE etaCIhi <- etaHat - etaCIlo etaCIlo <- etaHat + etaCIlo thetaHat <- etaHat thetaCIlo <- etaCIlo thetaCIhi <- etaCIhi methodParams <- NULL } else stop(paste("method",method,"not yet implemented")) } else if(fnctl=="geometric mean") { if(!is.null(replaceZeroes)) { u <- ya == 0 nReplace <- sum(u) if(is.logical(replaceZeroes)) replaceZeroes <- min(ya[!u]) ya[u] <- replaceZeroes } else if(any(ya <=0)) stop("cannot compute geometric mean with nonpositive data") if (method=="t.test") { mName <- "One sample t test on log transformed data" link <- "log" etaName <- fName etaHat <- mean(log(ya)) etaHatSE <- sqrt(var(log(ya)) / n) etaNull <- log(null.hypothesis) statistic <- (etaHat - etaNull) / etaHatSE df <- n - 1 if(test.type=="less") p.value <- pt(statistic,df) else if(test.type=="greater") p.value <- 1 - pt(statistic,df) else if(test.type=="two.sided") p.value <- 2 * pt(-abs(statistic),df) else p.value <- NA etaCIlo <- qt((1-conf.level)/2,df) * etaHatSE etaCIhi <- etaHat - etaCIlo etaCIlo <- etaHat + etaCIlo thetaHat <- exp(etaHat) thetaCIlo <- exp(etaCIlo) thetaCIhi <- exp(etaCIhi) methodParams <- NULL } else stop(paste("method",method,"not yet implemented")) } else if(fnctl=="proportion" || fnctl=="odds") { if (method=="KM") stop("method KM not yet implemented") else { if (method=="exact" || method=="exactLR") { binomInference <- binomInference.exactLR mName <- "exact distribution" if (test.type=="two.sided") mName <- paste(mName,"(LR ordering)") } else if (method=="exactTail") { binomInference <- binomInference.exactTail mName <- "exact distribution" if (test.type=="two.sided") mName <- paste(mName,"(tail probability ordering)") } else if (method=="wald") { binomInference <- binomInference.wald mName <- "Wald statistic" } else if (method=="cwald") { binomInference <- binomInference.cwald mName <- "continuity corrected Wald statistic" } else if (method=="score") { binomInference <- binomInference.score mName <- "score statistic" } else if (method=="cscore") { binomInference <- binomInference.cscore mName <- "continuity corrected score statistic" } else if (method=="agresti") { binomInference <- binomInference.agresti mName <- "Agresti & Coull" } else if (method=="jeffreys") { binomInference <- binomInference.jeffreys mName <- "Jeffreys" } else stop(paste("method",method,"not yet implemented")) mName <- paste("One sample inference for binomial proportions using", mName) link <- "identity" etaName <- fName z <- binomInference(y=sum(ya), n=length(ya), null.hypothesis=null.hypothesis, test.type=test.type, conf.level=conf.level) etaHat <- z[2] etaHatSE <- z[3] etaNull <- null.hypothesis statistic <- z[4] df <- NA p.value <- z[7] etaCIlo <- z[5] etaCIhi <- z[6] thetaHat <- etaHat thetaCIlo <- etaCIlo thetaCIhi <- etaCIhi methodParams <- NULL } } else stop(paste("inference for", fnctl, "not yet implemented")) } Inference <- cbind(n, thetaHat, thetaCIlo, thetaCIhi, null.hypothesis, p.value) dimnames(Inference) <- list("",c("n", fName, paste(format(100*conf.level),"% CIlo",sep=""), paste(format(100*conf.level),"% CIhi",sep=""), "Null Hyp", c("P", "P hi", "P lo", "P two")[tindx+1])) attr(Inference,"fnctl") <- fnctl attr(Inference,"hName") <- hName attr(Inference,"fName") <- fName attr(Inference,"method") <- method attr(Inference,"mName") <- mName attr(Inference,"methodParams") <- methodParams attr(Inference,"link") <- link attr(Inference,"isSurv") <- isSurv attr(Inference,"isDate") <- isDate attr(Inference,"isEvent") <- isEvent attr(Inference,"restriction") <- restriction attr(Inference,"rName") <- rName attr(Inference,"replaceZeroes") <- replaceZeroes attr(Inference,"nReplace") <- nReplace attr(Inference,"thresholds") <- thresholds attr(Inference,"nMsng") <- nMsng attr(Inference,"test.type") <- test.type attr(Inference,"conf.level") <- conf.level Statistics <- cbind(n, etaHat, etaHatSE, etaNull, statistic, df) dimnames(Statistics) <- list("",c("n", etaName, "SE", "Null", "TestStat", "df")) attr(Statistics,"link") <- link rslt <- list(Inference=Inference, Statistics=Statistics) class(rslt) <- "uOneSample" rslt } print.uOneSample <- function (obj,sigfigs=max(3,getOption("digits")-3),width=9,nonsci.limit=5, print.it= TRUE, version=F) { vrsn <- "20121107" if (version) return(vrsn) # # prints CI or exp(CI) to specified sigfigs # rounds exp(Est) to same number of decimal figures as exp(CI) # prints Std Err or Robust SE to specified sigfigs # rounds Estimate to same number of decimal figures as StdErr or Robust SE # prints F stat to specified sigfigs # prints P value to number of decimal figures as nonsci.limit unless < 10^-nonsci.limit when "< .00001" used # centers all but df and P value, which is right justified cmptRoundDigits <- function (x, sf) { y <- max(abs(x),na.rm=T) if (y==0) { sf } else { y <- trunc(log(y) / log(10)) - (y < 1) max(0,sf - y - (y < sf)) } } frmtCol <- function (x, sf, nonsci.limit, colwidth=9, append="") { rslt <- NULL for (i in 1:length(x)) { if (is.na(x[i])) { tmp <- "NA" } else { rd <- cmptRoundDigits (x[i], sf) if (rd <= nonsci.limit & abs(x[i]) < 10^nonsci.limit) { tmp <- format(round(x[i],rd),nsmall=rd,width=1) } else { tmp <- format(round(x[i],rd), digits=sf, scientific=T, width=1) } } rslt <- c(rslt,ifelse(x[i]<0,tmp,paste(" ",tmp,sep=""))) } rslt <- paste(rslt,append,sep="") format(rslt, justify="centre", width=colwidth) } x <- obj$Inference ncol <- dim(x)[2] frmtCoefficients <- format(x) frmtCoefficients[,1] <- format (x[,1],width=5) for (j in 2:ncol) frmtCoefficients[,j] <- frmtCol (x[,j],sigfigs,nonsci.limit,width) if (attr(x,"isDate") && attr(x,"fnctl") %in% c("mean","geometric mean","median","quantile")) { xformCol <- 2:4 orgn <- "1970-01-01" frmtCoefficients[x[,"isDate"]==1,xformCol] <- format(as.Date(x[x[,"isDate"]==1,xformCol],orgn)) } dimnames(frmtCoefficients)[[2]] <- format(dimnames(x)[[2]],justify="centre") restriction <- attr(x,"restriction") if (restriction < Inf) frmtCoefficients <- cbind(frmtCoefficients[,1,drop=F], "Restrict"=restriction,frmtCoefficients[,-1,drop=F]) if(print.it) { if (!is.null(attr(x,"hName"))) cat(attr(x,"hName"),"\n") cat("Method:",attr(x,"mName"),"\n") methodParams <- attr(x,"methodParams") if(!is.null(methodParams)) { cat(" Parameters:\n") print(methodParams) } nMsng <- attr(x, "nMsng") if(!is.na(nMsng) && nMsng > 0) cat(nMsng, "observations deleted due to missing values\n") print(frmtCoefficients,quote=F) } invisible(frmtCoefficients) }