Movatterモバイル変換

Version:	5.2-4
Date:	2025-10-02
Title:	Harrell Miscellaneous
Depends:	R (≥ 4.2.0)
Imports:	methods, ggplot2, cluster, rpart, nnet, foreign, gtable, grid,gridExtra, data.table, htmlTable (≥ 1.11.0), viridisLite,htmltools, base64enc, colorspace, rmarkdown, knitr, Formula
Suggests:	survival, qreport, acepack, chron, rms, mice, rstudioapi,tables, plotly (≥ 4.5.6), rlang, VGAM, leaps, pcaPP, digest,parallel, polspline, abind, kableExtra, rio, lattice,latticeExtra, gt, sparkline, jsonlite, htmlwidgets, qs,getPass, keyring, safer, htm2txt, boot
Description:	Contains many functions useful for dataanalysis, high-level graphics, utility operations, functions forcomputing sample size and power, simulation, importing and annotating datasets,imputing missing values, advanced table making, variable clustering,character string manipulation, conversion of R objects to LaTeX and html code,recoding variables, caching, simplified parallel computing, encrypting and decrypting data using a safe workflow, general moving window statistical estimation, and assistance in interpreting principal component analysis.
License:	GPL-2 |GPL-3 [expanded from: GPL (≥ 2)]
LazyLoad:	Yes
URL:	https://hbiostat.org/R/Hmisc/
Encoding:	UTF-8
RoxygenNote:	7.3.3
NeedsCompilation:	yes
Packaged:	2025-10-03 20:24:42 UTC; harrelfe
Author:	Frank E Harrell Jr [aut, cre], Cole Beck [ctb], Charles Dupont [ctb]
Maintainer:	Frank E Harrell Jr <fh@fharrell.com>
Repository:	CRAN
Date/Publication:	2025-10-05 06:50:02 UTC

Find Matching (or Non-Matching) Elements

Description

%nin% is a binary operator, which returns a logical vector indicatingif there is a match or not for its left operand. A true vector elementindicates no match in left operand, false indicates a match.

Usage

x %nin% table

Arguments

Value

vector of logical values with length equal to length ofx.

See Also

match%in%

Examples

c('a','b','c') %nin% c('a','b')

Character strings from unquoted names

Description

Cs makes a vector of character strings from a list of valid Rnames..q is similar but also makes uses of names of arguments.

Usage

Cs(...).q(...)

Arguments

Value

character string vector. For.q there will be anamesattribute to the vector if any names appeared in ....

See Also

sys.frame, deparse

Examples

Cs(a,cat,dog)# subset.data.frame <- dataframe[,Cs(age,sex,race,bloodpressure,height)].q(a, b, c, 'this and that').q(dog=a, giraffe=b, cat=c)

Empirical Cumulative Distribution Plot

Description

Computes coordinates of cumulative distribution function of x, and by defaultsplots it as a step function. A grouping variable may be specified so thatstratified estimates are computed and (by default) plotted. If there ismore than one group, thelabcurve function is used (by default) to labelthe multiple step functions or to draw a legend defining line types, colors,or symbols by linking them with group labels. Aweights vector maybe specified to get weighted estimates. Specifynormwt to makeweights sum to the length ofx (after removing NAs). Other wisethe total sample size is taken to be the sum of the weights.

Ecdf is actually a method, andEcdf.default is what'scalled for a vector argument.Ecdf.data.frame is called when thefirst argument is a data frame. This function can automatically set upa matrix of ECDFs and wait for a mouse click if the matrix requires morethan one page. Categorical variables, character variables, andvariables having fewer than a set number of unique values are ignored.Ifpar(mfrow=..) is not set up beforeEcdf.data.frame iscalled, the function will try to figure the best layout depending on thenumber of variables in the data frame. Upon return the originalmfrow is left intact.

When the first argument toEcdf is a formula, a Trellis/Lattice functionEcdf.formula is called. This allows for multi-panelconditioning, superposition using agroups variable, and otherTrellis features, along with the ability to easily plot transformedECDFs using thefun argument. For example, iffun=qnorm,the inverse normal transformation will be used for the y-axis. If thetransformed curves are linear this indicates normality. Like thexYplot function,Ecdf will create a functionKey ifthegroups variable is used. This function can be invoked by theuser to define the keys for the groups.

Usage

Ecdf(x, ...)## Default S3 method:Ecdf(x, what=c('F','1-F','f','1-f'),     weights=rep(1, length(x)), normwt=FALSE,     xlab, ylab, q, pl=TRUE, add=FALSE, lty=1,      col=1, group=rep(1,length(x)), label.curves=TRUE, xlim,      subtitles=TRUE, datadensity=c('none','rug','hist','density'),     side=1,      frac=switch(datadensity,none=NA,rug=.03,hist=.1,density=.1),     dens.opts=NULL, lwd=1, log='', ...)## S3 method for class 'data.frame'Ecdf(x, group=rep(1,nrows),     weights=rep(1, nrows), normwt=FALSE,     label.curves=TRUE, n.unique=10, na.big=FALSE, subtitles=TRUE,      vnames=c('labels','names'),...)## S3 method for class 'formula'Ecdf(x, data=sys.frame(sys.parent()), groups=NULL,     prepanel=prepanel.Ecdf, panel=panel.Ecdf, ..., xlab,     ylab, fun=function(x)x, what=c('F','1-F','f','1-f'), subset=TRUE)

Arguments

Value

forEcdf.default an invisible list with elements x and y giving thecoordinates of the cdf. If there is more than onegroup, a list ofsuch lists is returned. An attribute,N, is in the returnedobject. It contains the elementsn andm, the number ofnon-missing and missing observations, respectively.

Side Effects

plots

Author(s)

Frank Harrell
Department of Biostatistics, Vanderbilt University
fh@fharrell.com

See Also

wtd.Ecdf,label,table,cumsum,labcurve,xYplot,histSpike

Examples

set.seed(1)ch <- rnorm(1000, 200, 40)Ecdf(ch, xlab="Serum Cholesterol")scat1d(ch)                       # add rug plothistSpike(ch, add=TRUE, frac=.15)   # add spike histogram# Better: add a data density display automatically:Ecdf(ch, datadensity='density')label(ch) <- "Serum Cholesterol"Ecdf(ch)other.ch <- rnorm(500, 220, 20)Ecdf(other.ch,add=TRUE,lty=2)sex <- factor(sample(c('female','male'), 1000, TRUE))Ecdf(ch, q=c(.25,.5,.75))  # show quartilesEcdf(ch, group=sex,     label.curves=list(method='arrow'))# Example showing how to draw multiple ECDFs from paired datapre.test <- rnorm(100,50,10)post.test <- rnorm(100,55,10)x <- c(pre.test, post.test)g <- c(rep('Pre',length(pre.test)),rep('Post',length(post.test)))Ecdf(x, group=g, xlab='Test Results', label.curves=list(keys=1:2))# keys=1:2 causes symbols to be drawn periodically on top of curves# Draw a matrix of ECDFs for a data framem <- data.frame(pre.test, post.test,                 sex=sample(c('male','female'),100,TRUE))Ecdf(m, group=m$sex, datadensity='rug')freqs <- sample(1:10, 1000, TRUE)Ecdf(ch, weights=freqs)  # weighted estimates# Trellis/Lattice examples:region <- factor(sample(c('Europe','USA','Australia'),100,TRUE))year <- factor(sample(2001:2002,1000,TRUE))Ecdf(~ch | region*year, groups=sex)Key()           # draw a key for sex at the default location# Key(locator(1)) # user-specified positioning of keyage <- rnorm(1000, 50, 10)Ecdf(~ch | lattice::equal.count(age), groups=sex)  # use overlapping shinglesEcdf(~ch | sex, datadensity='hist', side=3)  # add spike histogram at top

Debug Printing Function Generator

Description

Takes the name of a systemoptions(opt=) and checks to see if optionopt isset toTRUE, taking its default value to beFALSE. IfTRUE, a function iscreated that callsprn() to print an object with the object's name in thedescription along with the option name and the name of the function within whichthe generated function was called, if any. If optionopt is not set, a dummy functionis generated instead. Ifoptions(debug_file=) is set when the generated functionis called,prn() output will be appended to that file name instead of the console.At any time, setoptions(debug_file='') to resume printing to the console.

Usage

Fdebug(opt)

Arguments

Value

a function

Author(s)

Fran Harrell

Examples

dfun <- Fdebug('my_option_name')   # my_option_name not currently setdfundfun(sqrt(2))options(my_option_name=TRUE)dfun <- Fdebug('my_option_name')dfundfun(sqrt(2))# options(debug_file='/tmp/z') to append output to /tmp/zoptions(my_option_name=NULL)

Gini's Mean Difference

Description

GiniMD computes Gini's mean difference on anumeric vector. This index is defined as the mean absolute differencebetween any two distinct elements of a vector. For a Bernoulli(binary) variable with proportion of ones equal top and samplesizen, Gini's mean difference is2\frac{n}{n-1}p(1-p). For a trinomial variable (e.g., predicted values for a 3-level categoricalpredictor using two dummy variables) having (predicted)valuesA, B, C with corresponding proportionsa, b, c,Gini's mean difference is2\frac{n}{n-1}[ab|A-B|+ac|A-C|+bc|B-C|]

Usage

GiniMd(x, na.rm=FALSE)

Arguments

Value

a scalar numeric

Author(s)

Frank Harrell
Department of Biostatistics
Vanderbilt University
fh@fharrell.com

References

David HA (1968): Gini's mean difference rediscovered. Biometrika 55:573–575.

Examples

set.seed(1)x <- rnorm(40)# Test GiniMd against a brute-force solutiongmd <- function(x) {  n <- length(x)  sum(outer(x, x, function(a, b) abs(a - b))) / n / (n - 1)  }GiniMd(x)gmd(x)z <- c(rep(0,17), rep(1,6))n <- length(z)GiniMd(z)2*mean(z)*(1-mean(z))*n/(n-1)a <- 12; b <- 13; c <- 7; n <- a + b + cA <- -.123; B <- -.707; C <- 0.523xx <- c(rep(A, a), rep(B, b), rep(C, c))GiniMd(xx)2*(a*b*abs(A-B) + a*c*abs(A-C) + b*c*abs(B-C))/n/(n-1)

Internal Hmisc functions

Description

Internal Hmisc functions.

Details

These are not to be called by the user or are undocumented.

Overview of Hmisc Library

Description

The Hmisc library contains many functions useful for data analysis,high-level graphics, utility operations, functions for computingsample size and power, translating SAS datasets intoR, imputingmissing values, advanced table making, variable clustering, characterstring manipulation, conversion ofR objects to LaTeX code, recodingvariables, and bootstrap repeated measures analysis. Most of thesefunctions were written by F Harrell, but a few were collected fromstatlib and from s-news; other authors are indicated below. Thiscollection of functions includes all of Harrell's submissions tostatlib other than the functions in therms and displaylibraries. A few of the functions do not have “Help”documentation.

To makeHmisc load silently, issueoptions(Hverbose=FALSE) beforelibrary(Hmisc).

Functions

Copyright Notice

GENERAL DISCLAIMER
This program is free software; you can redistribute itand/or modify it under the terms of the GNU General PublicLicense as published by the Free Software Foundation; eitherversion 2, or (at your option) any later version.

This program is distributed in the hope that it will beuseful, but WITHOUT ANY WARRANTY; without even the impliedwarranty of MERCHANTABILITY or FITNESS FOR A PARTICULARPURPOSE. See the GNU General Public License for moredetails.

In short: You may use it any way you like, as long as youdon't charge money for it, remove this notice, or hold anyone liablefor its results. Also, please acknowledge the source and communicatechanges to the author.

If this software is used is work presented for publication, kindlyreference it using for example:
Harrell FE (2014): Hmisc: A package of miscellaneous R functions.Programs available fromhttps://hbiostat.org/R/Hmisc/.
Be sure to referenceR itself and other libraries used.

Author(s)

Frank E Harrell Jr
Professor of Biostatistics
Vanderbilt University School of Medicine
Nashville, Tennessee
fh@fharrell.com

References

See Alzola CF, Harrell FE (2004): An Introduction to S and theHmisc and Design Libraries athttps://hbiostat.org/R/doc/sintro.pdffor extensive documentation and examples for the Hmisc package.

Lag a Numeric, Character, or Factor Vector

Description

Shifts a vectorshift elements later. Character or factorvariables are padded with"", numerics withNA. The shiftmay be negative.

Usage

Lag(x, shift = 1)

Arguments

Details

A.ttributes of the original object are carried along to the new laggedone.

Value

a vector likex

Author(s)

Frank Harrell

See Also

lag

Examples

Lag(1:5,2)Lag(letters[1:4],2)Lag(factor(letters[1:4]),-2)# Find which observations are the first for a given subjectid <- c('a','a','b','b','b','c')id != Lag(id)!duplicated(id)

Merge Multiple Data Frames or Data Tables

Description

Merges an arbitrarily large series of data frames or data tables containing commonid variables. Information about number of observations and number of uniqueids in individual and final merged datasets is printed. The first data frame/table has special meaning in that all of its observations are kept whether they matchids in other data frames or not. For all other data frames, by default non-matching observations are dropped. The first data frame is also the one against which counts of uniqueids are compared. Sometimesmerge drops variable attributes such aslabels andunits. These are restored byMerge.

Usage

Merge(..., id = NULL, all = TRUE, verbose = TRUE)

Arguments

Examples

## Not run: a <- data.frame(sid=1:3, age=c(20,30,40))b <- data.frame(sid=c(1,2,2), bp=c(120,130,140))d <- data.frame(sid=c(1,3,4), wt=c(170,180,190))all <- Merge(a, b, d, id = ~ sid)# First file should be the master file and must# contain all ids that ever occur.  ids not in the master will# not be merged from other datasets.a <- data.table(a); setkey(a, sid)# data.table also does not allow duplicates without allow.cartesian=TRUEb <- data.table(sid=1:2, bp=c(120,130)); setkey(b, sid)d <- data.table(d); setkey(d, sid)all <- Merge(a, b, d)## End(Not run)

Miscellaneous Functions

Description

This documents miscellaneous small functions in Hmisc that may be ofinterest to users.

clowess runslowess but if theiter argumentexceeds zero, sometimes wild values can result, in which caselowess is re-run withiter=0.

confbar draws multi-level confidence bars using small rectanglesthat may be of different colors.

getLatestSource fetches andsources the most recentsource code for functions in GitHub.

grType retrieves the system optiongrType, which isforced to be"base" if theplotly package is notinstalled.

prType retrieves the system optionprType, which isset to"plain" if the option is not set.print methodsthat allow for markdown/html/latex can be automatically invoked bysettingoptions(prType="html") oroptions(prType='latex').

htmlSpecialType retrieves the system optionhtmlSpecialType, which is set to"unicode" if the optionis not set.htmlSpecialType='unicode' cause html-generatingfunctions inHmisc andrms to use unicode for specialcharacters, andhtmlSpecialType='&' uses the older ampersand3-digit format.

inverseFunction generates a function to find all inverses of amonotonic or nonmonotonic function that is tabulated at vectors (x,y),typically 1000 points. If the original function is monotonic, simple linearinterpolation is used and the result is a vector, otherwise linearinterpolation is used within each interval in which the function ismonotonic and the result is a matrix with number of columns equal to thenumber of monotonic intervals. If a requested y is not within anyinterval, the extreme x that pertains to the nearest extreme y isreturned. Specifying what='sample' to the returned function will cause avector to be returned instead of a matrix, with elements taken as arandom choice of the possible inverses.

james.stein computes James-Stein shrunken estimates of cellmeans given a response variable (which may be binary) and a groupingindicator.

keepHattrib for an input variable or a data frame, creates alist object saving special Hmisc attributes such aslabel andunits that might be lost during certain operations such asrunningdata.table.restoreHattrib restores these attributes.

km.quick provides a fast way to invokesurvfitKM in thesurvival package to efficiently get Kaplan-Meier or Fleming-Harrington estimates for asingle stratum for a vector of time points (iftimes is given) or toget a vector of survival time quantiles (ifq is given). If neither is given,the whole curve is returned in a list with objectstime andsurv, andthere is an option to consider an interval as pertaining to greater than or equalto a specific time instead of the traditional greater than. If the censoring is not right censoring, the more generalsurvfit is called bykm.quick.

latexBuild takes pairs of character strings and produces asingle character string containing concatenation of all of them, plusan attribute"close" which is a character string containing theLaTeX closure that will balance LaTeX code with respect toparentheses, braces, brackets, orbegin vs.end. Whenan even-numbered element of the vector is not a left parenthesis,brace, or bracket, the element is taken as a word that was surroundedbybegin and braces, for which the correspondingend isconstructed in the returned attribute.

lm.fit.qr.bare is a fast stripped-down function for computingregression coefficients, residuals,R^2, and fitted values. Ituseslm.fit.

matxv multiplies a matrix by a vector, handling automaticaddition of intercepts if the matrix does not have a column of ones.If the first argument is not a matrix, it will be converted to one.An optional argument allows the second argument to be treated as amatrix, useful when its rows represent bootstrap reps ofcoefficients. Then ab' is computed.matxv respects the"intercepts" attribute if it is stored onb by therms package. This is used byormfits that are bootstrap-repeated bybootcov whereonly the intercept corresponding to the median is retained. Ifkint has nonzero length, it is checked for consistency with theattribute.

makeSteps is a copy of the dostep function inside thesurvival package'splot.survfit function. It expands aseries of points to include all the segments needed to plot stepfunctions. This is useful for drawing polygons to shade confidencebands for step functions.

nomiss returns a data frame (if its argument is one) with rowscorresponding toNAs removed, or it returns a matrix with rowswith any element missing removed.

outerText usesaxis() to put right-justified textstrings in the right margin. Placement depends onpar('mar')[4]

plotlyParm is a list of functions useful for specifyingparameters toplotly graphics.

plotp is a generic to handleplotp methods to makeplotly graphics.

rendHTML renders HTML in a character vector, first convertingto one character string with newline delimeters. Ifknitr iscurrently running, runs this string throughknitr::asis_outputso that the user need not includeresults='asis' in the chunkheader for R Markdown or Quarto. Ifknitr is not running, useshtmltools::browsable andhtmltools::HTML and prints theresult so that an RStudio viewer (if running inside RStudio) orseparate browser window displays the rendered HTML. The HTML code issurrounded by yaml markup to make Pandoc not fiddle with the HTML.Set the argumenthtml=FALSE to not add this, in case you arereally rendering markdown.html=FALSE also invokesrmarkdown::render to convert the character vector to HTMLbefore usinghtmltools to view, assuming the charactersrepresent RMarkdown/Quarto text other than the YAML header. Ifoptions(rawmarkup=TRUE) is in effect,rendHTML will justcat() its first argument. This is useful when rendering ishappening inside a Quarto margin, for example.

sepUnitsTrans converts character vectors containing values suchasc("3 days","3day","4month","2 years","2weeks","7") tonumeric vectors (herec(3,3,122,730,14,7)) in a flexible fashion. The user canspecify a vector of units of measurements and conversion factors. The unitswith a conversion factor of1 are taken as the target units,and if those units are present in the character strings they areignored. The target units are added to the resulting vector as the"units" attribute.

strgraphwrap is likestrwrap but is for the currentgraphics environment.

tobase64image is a function written by Dirk Eddelbuettel thatuses thebase64enc package to convert a png graphic file tobase64 encoding to include as an inline image in an html file.

trap.rule computes the area under a curve using the trapezoidalrule, assumingx is sorted.

trellis.strip.blank sets up Trellis or Lattice graphs to have aclear background on the strips for panel labels.

unPaste provides a version of the S-Plusunpaste thatworks forR and S-Plus.

whichClosePW is a very fast function using weighted multinomialsampling to determine which element of a vector is "closest" to eachelement of another vector.whichClosest quickly finds the closestelement without any randomness.

whichClosek is a slow function that finds, after jittering thelookup table, thek closest matchest to each element of theother vector, and chooses from among these one at random.

xless is a function for Linux/Unix users to invoke the systemxless command to pop up a window to display the result ofprinting an object. For MacOSxless uses the systemopen command to pop up aTextEdit window.

Usage

confbar(at, est, se, width, q = c(0.7, 0.8, 0.9, 0.95, 0.99),         col = gray(c(0, 0.25, 0.5, 0.75, 1)),        type = c("v", "h"), labels = TRUE, ticks = FALSE,        cex = 0.5, side = "l", lwd = 5, clip = c(-1e+30, 1e+30),        fun = function(x) x,        qfun = function(x) ifelse(x == 0.5, qnorm(x),                            ifelse(x < 0.5, qnorm(x/2),                            qnorm((1 +  x)/2))))getLatestSource(x=NULL, package='Hmisc', recent=NULL, avail=FALSE)grType()prType()htmlSpecialType()inverseFunction(x, y)james.stein(y, group)keepHattrib(obj)km.quick(S, times, q,        type = c("kaplan-meier", "fleming-harrington", "fh2"),        interval = c(">", ">="), method=c('constant', 'linear'), fapprox=0, n.risk=FALSE)latexBuild(..., insert, sep='')lm.fit.qr.bare(x, y, tolerance, intercept=TRUE, xpxi=FALSE, singzero=FALSE)matxv(a, b, kint=1, bmat=FALSE)nomiss(x)outerText(string, y, cex=par('cex'), ...)plotlyParmplotp(data, ...)rendHTML(x, html=TRUE)restoreHattrib(obj, attribs)sepUnitsTrans(x, conversion=c(day=1, month=365.25/12, year=365.25, week=7),              round=FALSE, digits=0)strgraphwrap(x, width = 0.9 * getOption("width"),             indent = 0, exdent = 0,             prefix = "", simplify = TRUE, units='user', cex=NULL)tobase64image(file, Rd = FALSE, alt = "image")trap.rule(x, y)trellis.strip.blank()unPaste(str, sep="/")whichClosest(x, w)whichClosePW(x, w, f=0.2)whichClosek(x, w, k)xless(x, ..., title)

Arguments

Author(s)

Frank Harrell and Charles Dupont

Examples

trap.rule(1:100,1:100)unPaste(c('a;b or c','ab;d','qr;s'), ';')sepUnitsTrans(c('3 days','4 months','2 years','7'))set.seed(1)whichClosest(1:100, 3:5)whichClosest(1:100, rep(3,20))whichClosePW(1:100, rep(3,20))whichClosePW(1:100, rep(3,20), f=.05)whichClosePW(1:100, rep(3,20), f=1e-10)x <- seq(-1, 1, by=.01)y <- x^2h <- inverseFunction(x,y)formals(h)$turns   # vertexa <- seq(0, 1, by=.01)plot(0, 0, type='n', xlim=c(-.5,1.5))lines(a, h(a)[,1])            ## first inverselines(a, h(a)[,2], col='red') ## second inversea <- c(-.1, 1.01, 1.1, 1.2)points(a, h(a)[,1])d <- data.frame(x=1:2, y=3:4, z=5:6)d <- upData(d, labels=c(x='X', z='Z lab'), units=c(z='mm'))a <- keepHattrib(d)d <- data.frame(x=1:2, y=3:4, z=5:6)d2 <- restoreHattrib(d, a)sapply(d2, attributes)## Not run: getLatestSource(recent=5)  # source() most recent 5 revised files in HmiscgetLatestSource('cut2')    # fetch and source latest cut2.sgetLatestSource('all')     # get everythinggetLatestSource(avail=TRUE) # list available files and latest versions## End(Not run)

R2Measures

Description

Generalized R^2 Measures

Usage

R2Measures(lr, p, n, ess = NULL, padj = 1)

Arguments

Details

Computes various generalized R^2 measures related to the Maddala-Cox-Snell (MCS) R^2 for regression models fitted with maximum likelihood. The original MCS R^2 is labeledR2 in the result. This measure uses the raw sample sizen and does not penalize for the number of free parameters, so it can be rewarded for overfitting. A measure adjusted for the number of fitted regression coefficientsp uses the analogy to R^2 in linear models by computing1 - exp(- lr / n) * (n-1)/(n-p-1) ifpadj=2, which is approximately1 - exp(- (lr - p) / n), the version used ifpadj=1 (the default). The latter measure is appealing because the expected value of the likelihood ratio chi-square statisticlr isp under the global null hypothesis of no predictors being associated with the response variable. Seehttps://hbiostat.org/bib/r2.html for more details.

It is well known that in logistic regression the MCS R^2 cannot achieve a value of 1.0 even with a perfect model, which prompted Nagelkerke to divide the R^2 measure by its maximum attainable value. This is not necessarily the best recalibration of R^2 throughout its range. An alternative is to use the formulas above but to replace the raw sample sizen with the effective sample size, which for data with many ties can be significantly lower than the number of observations. As used in thepopower() anddescribe() functions, in the context of a Wilcoxon test or the proportional odds model, the effective sample size isn * (1 - f) wheref is the sums of cubes of the proportion of observations at each distict value of the response variable. Whitehead derived this from an approximation to the variance of a log odds ratio in a proportional odds model. To obtain R^2 measures using the effective sample size, either provideess as a single number specifying the effective sample size, or specify a vector of frequencies of distinct Y values from which the effective sample size will be computed. In the context of survival analysis, the single number effective sample size you may wish to specify is the number of uncensored observations. This is exactly correct when estimating the hazard rate from a simple exponential distribution or when using the Cox PH/log-rank test. For failure time distributions with a very high early hazard, censored observations contain enough information that the effective sample size is greater than the number of events. See Benedetti et al, 1982.

If the effective sample size equals the raw sample size, measures involving the effective sample size are set toNA.

Value

named vector of R2 measures. The notation for results isR^2(p, n) where thep component is empty for unadjusted estimates andn is the sample size used (actual sample size for first measures, effective sample size for remaining ones). For indexes that are not adjusted, onlyn appears.

Author(s)

Frank Harrell

References

Smith TJ and McKenna CM (2013): A comparison of logistic regression pseudo R^2 indices. Multiple Linear Regression Viewpoints 39:17-26.https://www.glmj.org/archives/articles/Smith_v39n2.pdf

Benedetti JK, et al (1982): Effective sample size for tests of censored survival data. Biometrika 69:343–349.

Mittlbock M, Schemper M (1996): Explained variation for logistic regression. Stat in Med 15:1987-1997.

Date, S: R-squared, adjusted R-squared and pseudo R-squared.https://timeseriesreasoning.com/contents/r-squared-adjusted-r-squared-pseudo-r-squared/

UCLA: What are pseudo R-squareds?https://stats.oarc.ucla.edu/other/mult-pkg/faq/general/faq-what-are-pseudo-r-squareds/

Allison P (2013): What's the beset R-squared for logistic regression?https://statisticalhorizons.com/r2logistic/

Menard S (2000): Coefficients of determination for multiple logistic regression analysis. The Am Statistician 54:17-24.

Whitehead J (1993): Sample size calculations for ordered categorical data. Stat in Med 12:2257-2271. See errata (1994) 13:871 and letter to the editor by Julious SA, Campbell MJ (1996) 15:1065-1066 showing that for 2-category Y the Whitehead sample size formula agrees closely with the usual formula for comparing two proportions.

Examples

x <- c(rep(0, 50), rep(1, 50))y <- x# f <- lrm(y ~ x)# f   # Nagelkerke R^2=1.0# lr <- f$stats['Model L.R.']# 1 - exp(- lr / 100)  # Maddala-Cox-Snell (MCS) 0.75lr <- 138.6267  # manually so don't need rms packageR2Measures(lr, 1, 100, c(50, 50))  # 0.84 Effective n=75R2Measures(lr, 1, 100, 50)         # 0.94# MCS requires unreasonable effective sample size = minimum outcome# frequency to get close to the 1.0 that Nagelkerke R^2 achieves

Faciliate Use of save and load to Remote Directories

Description

These functions are slightly enhanced versions ofsave andload that allow a target directory to be specified usingoptions(LoadPath="pathname"). If theLoadPath option isnot set, the current working directory is used.

Usage

# options(LoadPath='mypath')Save(object, name=deparse(substitute(object)), compress=TRUE)Load(object)

Arguments

Details

Save creates a temporary version of the object under the namegiven by the user, so thatsave will internalize this name.Then subsequentLoad orload will cause an object of theoriginal name to be created in the global environment. The name oftheR data file is assumed to be the name of the object (or the valueofname) appended with".rda".

Author(s)

Frank Harrell

See Also

save,load

Examples

## Not run: d <- data.frame(x=1:3, y=11:13)options(LoadPath='../data/rda')Save(d)   # creates ../data/rda/d.rdaLoad(d)   # reads   ../data/rda/d.rdaSave(d, 'D')   # creates object D and saves it in .../D.rda## End(Not run)

Indexes of Absolute Prediction Error for Linear Models

Description

Computes the mean and median of various absolute errors related toordinary multiple regression models. The mean and median absoluteerrors correspond to the mean square due to regression, error, andtotal. The absolute errors computed are derived from\hat{Y} - \mbox{median($\hat{Y}$)},\hat{Y} - Y, andY - \mbox{median($Y$)}. The function alsocomputes ratios that correspond toR^2 and1 - R^2 (butthese ratios do not add to 1.0); theR^2 measure is the ratio ofmean or median absolute\hat{Y} - \mbox{median($\hat{Y}$)} to the mean or median absoluteY - \mbox{median($Y$)}. The1 - R^2 or SSE/SSTmeasure is the mean or median absolute\hat{Y} - Ydivided by the mean or median absolute\hat{Y} - \mbox{median($Y$)}.

Usage

abs.error.pred(fit, lp=NULL, y=NULL)## S3 method for class 'abs.error.pred'print(x, ...)

Arguments

Value

a list of classabs.error.pred (used byprint.abs.error.pred) containing two matrices:differences andratios.

Author(s)

Frank Harrell
Department of Biostatistics
Vanderbilt University School of Medicine
fh@fharrell.com

References

Schemper M (2003): Stat in Med 22:2299-2308.

Tian L, Cai T, Goetghebeur E, Wei LJ (2007): Biometrika 94:297-311.

See Also

lm,ols,cor,validate.ols

Examples

set.seed(1)         # so can regenerate resultsx1 <- rnorm(100)x2 <- rnorm(100)y  <- exp(x1+x2+rnorm(100))f <- lm(log(y) ~ x1 + poly(x2,3), y=TRUE)abs.error.pred(lp=exp(fitted(f)), y=y)rm(x1,x2,y,f)

Add Marginal Observations

Description

Given a data frame and the names of variable, doubles thedata frame for each variable with a new category"All" by default, or by the value oflabel.A new variable.marginal. is added to the resulting data frame,with value"" if the observation is an original one, and withvalue equal to the names of the variable being marginalized (separatedby commas) otherwise. If there is another stratification variablebesides the one in ..., and that variable is nested inside thevariable in ..., specifynested=variable name to have the valueof that variable set folabel whenever marginal observations arecreated for .... See the state-city example below.

Usage

addMarginal(data, ..., label = "All", margloc=c('last', 'first'), nested)

Arguments

Examples

d <- expand.grid(sex=c('female', 'male'), country=c('US', 'Romania'),                 reps=1:2)addMarginal(d, sex, country)# Example of nested variablesd <- data.frame(state=c('AL', 'AL', 'GA', 'GA', 'GA'),                city=c('Mobile', 'Montgomery', 'Valdosto',                       'Augusta', 'Atlanta'),                x=1:5, stringsAsFactors=TRUE)addMarginal(d, state, nested=city) # cite set to 'All' when state is

addggLayers

Description

Add Spike Histograms and Extended Box Plots toggplot

Usage

addggLayers(  g,  data,  type = c("ebp", "spike"),  ylim = layer_scales(g)$y$get_limits(),  by = "variable",  value = "value",  frac = 0.065,  mult = 1,  facet = NULL,  pos = c("bottom", "top"),  showN = TRUE)

Arguments

Details

For an example seethis. Note that it was not possible to just create the layers needed to be added, as creating these particular layers in isolation resulted in aggplot error.

Value

the originalggplot object with more layers added

Author(s)

Frank Harrell

See Also

spikecomp()

Check if All Elements in Character Vector are Numeric

Description

Tests, without issuing warnings, whether all elements of a charactervector are legal numeric values, or optionally converts the vector to anumeric vector. Leading and trailing blanks inx are ignored.

Usage

all.is.numeric(x, what = c("test", "vector", "nonnum"), extras=c('.','NA'))

Arguments

Value

a logical value ifwhat="test" or a vector otherwise

Author(s)

Frank Harrell

See Also

as.numeric

Examples

all.is.numeric(c('1','1.2','3'))all.is.numeric(c('1','1.2','3a'))all.is.numeric(c('1','1.2','3'),'vector')all.is.numeric(c('1','1.2','3a'),'vector')all.is.numeric(c('1','',' .'),'vector')all.is.numeric(c('1', '1.2', '3a'), 'nonnum')

Linear Extrapolation

Description

Works in conjunction with theapprox function to do linearextrapolation.approx in R does not support extrapolation atall, and it is buggy in S-Plus 6.

Usage

approxExtrap(x, y, xout, method = "linear", n = 50, rule = 2, f = 0,             ties = "ordered", na.rm = FALSE)

Arguments

Details

Duplicates inx (and correspondingy elements) are removedbefore usingapprox.

Value

a vector the same length asxout

Author(s)

Frank Harrell

See Also

approx

Examples

approxExtrap(1:3,1:3,xout=c(0,4))

Additive Regression with Optimal Transformations on Both Sides usingCanonical Variates

Description

Expands continuous variables into restricted cubic spline bases andcategorical variables into dummy variables and fits a multivariateequation using canonical variates. This finds optimum transformationsthat maximizeR^2. Optionally, the bootstrap is used to estimatethe covariance matrix of both left- and right-hand-side transformationparameters, and to estimate the bias in theR^2 due to overfittingand compute the bootstrap optimism-correctedR^2.Cross-validation can also be used to get an unbiased estimate ofR^2 but this is not as precise as the bootstrap estimate. Thebootstrap and cross-validation may also used to get estimates of meanand median absolute error in predicted values on the originalyscale. These two estimates are perhaps the best ones for gauging theaccuracy of a flexible model, because it is difficult to compareR^2 under different y-transformations, and becauseR^2allows for an out-of-sample recalibration (i.e., it only measuresrelative errors).

Note that uncertainty about the proper transformation ofy causesan enormous amount of model uncertainty. When the transformation fory is estimated from the data a high variance in predicted valueson the originaly scale may result, especially if the truetransformation is linear. Comparing bootstrap or cross-validated meanabsolute errors with and without restricted they transform to belinear (ytype='l') may help the analyst choose the proper modelcomplexity.

Usage

areg(x, y, xtype = NULL, ytype = NULL, nk = 4,     B = 0, na.rm = TRUE, tolerance = NULL, crossval = NULL)## S3 method for class 'areg'print(x, digits=4, ...)## S3 method for class 'areg'plot(x, whichx = 1:ncol(x$x), ...)## S3 method for class 'areg'predict(object, x, type=c('lp','fitted','x'),                       what=c('all','sample'), ...)

Arguments

Details

areg is a competitor oface in theacepackpackage. Transformations fromace are seldom smooth enough andare often overfitted. Withareg the complexity can be controlledwith thenk parameter, and predicted values are easy to obtainbecause parametric functions are fitted.

If one side of the equation has a categorical variable with more thantwo categories and the other side has a continuous variable not assumedto act linearly, larger sample sizes are needed to reliably estimatetransformations, as it is difficult to optimally score categoricalvariables to maximizeR^2 against a simultaneously optimallytransformed continuous variable.

Value

a list of class"areg" containing many objects

Author(s)

Frank Harrell
Department of Biostatistics
Vanderbilt University
fh@fharrell.com

References

Breiman and Friedman, Journal of the American StatisticalAssociation (September, 1985).

See Also

cancor,ace,transcan

Examples

set.seed(1)ns <- c(30,300,3000)for(n in ns) {  y <- sample(1:5, n, TRUE)  x <- abs(y-3) + runif(n)  par(mfrow=c(3,4))  for(k in c(0,3:5)) {    z <- areg(x, y, ytype='c', nk=k)    plot(x, z$tx)title(paste('R2=',format(z$rsquared)))    tapply(z$ty, y, range)    a <- tapply(x,y,mean)    b <- tapply(z$ty,y,mean)    plot(a,b)abline(lsfit(a,b))    # Should get same result to within linear transformation if reverse x and y    w <- areg(y, x, xtype='c', nk=k)    plot(z$ty, w$tx)    title(paste('R2=',format(w$rsquared)))    abline(lsfit(z$ty, w$tx)) }}par(mfrow=c(2,2))# Example where one category in y differs from others but only in variance of xn <- 50y <- sample(1:5,n,TRUE)x <- rnorm(n)x[y==1] <- rnorm(sum(y==1), 0, 5)z <- areg(x,y,xtype='l',ytype='c')zplot(z)z <- areg(x,y,ytype='c')zplot(z)## Not run: # Examine overfitting when true transformations are linearpar(mfrow=c(4,3))for(n in c(200,2000)) {  x <- rnorm(n); y <- rnorm(n) + x    for(nk in c(0,3,5)) {    z <- areg(x, y, nk=nk, crossval=10, B=100)    print(z)    plot(z)    title(paste('n=',n))  }}par(mfrow=c(1,1))# Underfitting when true transformation is quadratic but overfitting# when y is allowed to be transformedset.seed(49)n <- 200x <- rnorm(n); y <- rnorm(n) + .5*x^2#areg(x, y, nk=0, crossval=10, B=100)#areg(x, y, nk=4, ytype='l', crossval=10, B=100)z <- areg(x, y, nk=4) #, crossval=10, B=100)z# Plot x vs. predicted value on original scale.  Since y-transform is# not monotonic, there are multiple y-inversesxx <- seq(-3.5,3.5,length=1000)yhat <- predict(z, xx, type='fitted')plot(x, y, xlim=c(-3.5,3.5))for(j in 1:ncol(yhat)) lines(xx, yhat[,j], col=j)# Plot a random sample of possible y inversesyhats <- predict(z, xx, type='fitted', what='sample')points(xx, yhats, pch=2)## End(Not run)# True transformation of x1 is quadratic, y is linearn <- 200x1 <- rnorm(n); x2 <- rnorm(n); y <- rnorm(n) + x1^2z <- areg(cbind(x1,x2),y,xtype=c('s','l'),nk=3)par(mfrow=c(2,2))plot(z)# y transformation is inverse quadratic but areg gets the same answer by# making x1 quadraticn <- 5000x1 <- rnorm(n); x2 <- rnorm(n); y <- (x1 + rnorm(n))^2z <- areg(cbind(x1,x2),y,nk=5)par(mfrow=c(2,2))plot(z)# Overfit 20 predictors when no true relationships existn <- 1000x <- matrix(runif(n*20),n,20)y <- rnorm(n)z <- areg(x, y, nk=5)  # add crossval=4 to expose the problem# Test predict functionn <- 50x <- rnorm(n)y <- rnorm(n) + xg <- sample(1:3, n, TRUE)z <- areg(cbind(x,g),y,xtype=c('s','c'))range(predict(z, cbind(x,g)) - z$linear.predictors)

Multiple Imputation using Additive Regression, Bootstrapping, andPredictive Mean Matching

Description

Thetranscan function creates flexible additive imputation modelsbut provides only an approximation to true multiple imputation as theimputation models are fixed before all multiple imputations aredrawn. This ignores variability caused by having to fit theimputation models.aregImpute takes all aspects of uncertainty inthe imputations into account by using the bootstrap to approximate theprocess of drawing predicted values from a full Bayesian predictivedistribution. Different bootstrap resamples are used for each of themultiple imputations, i.e., for theith imputation of a sometimesmissing variable,i=1,2,... n.impute, a flexible additivemodel is fitted on a sample with replacement from the original data andthis model is used to predict all of the original missing andnon-missing values for the target variable.

areg is used to fit the imputation models. By default, linearityis assumed for target variables (variables being imputed) andnk=3 knots are assumed for continuous predictors transformedusing restricted cubic splines. Ifnk is three or greater andtlinear is set toFALSE,aregsimultaneously finds transformations of the target variable and of all ofthe predictors, to get a good fit assuming additivity, maximizingR^2, using the same canonical correlation method astranscan. Flexible transformations may be overridden forspecific variables by specifying the identity transformation for them.When a categorical variable is being predicted, the flexibletransformation is Fisher's optimum scoring method. Nonlinear transformations for continuous variables may be nonmonotonic. Ifnk is a vector,areg's bootstrap andcrossval=10options will be used to help find the optimum validating value ofnk over values of that vector, at the last imputation iteration.For the imputations, the minimum value ofnk is used.

Instead of defaulting to taking random draws from fitted imputationmodels using random residuals as is done bytranscan,aregImpute by default uses predictive mean matching with optionalweighted probability sampling of donors rather than using only theclosest match. Predictive mean matching works for binary, categorical,and continuous variables without the need for iterative maximumlikelihood fitting for binary and categorical variables, and without theneed for computing residuals or for curtailing imputed values to be inthe range of actual data. Predictive mean matching is especiallyattractive when the variable being imputed is also being transformedautomatically. Constraints may be placed on variables being imputedwith predictive mean matching, e.g., a missing hospital discharge datemay be required to be imputed from a donor observation whose dischargedate is before the recipient subject's first post-discharge visit date.See Details below for more information about thealgorithm. A"regression" method is also available that issimilar to that used intranscan. This option should be usedwhen mechanistic missingness requires the use of extrapolation duringimputation.

Aprint method summarizes the results, and aplot method plotsdistributions of imputed values. Typically,fit.mult.impute willbe called afteraregImpute.

If a target variable is transformed nonlinearly (i.e., ifnk isgreater than zero andtlinear is set toFALSE) and theestimated target variable transformation is non-monotonic, imputedvalues are not unique. Whentype='regression', a random choiceof possible inverse values is made.

ThereformM function provides two ways of recreating a formula togive toaregImpute by reordering the variables in the formula.This is a modified version of a function written by Yong Hao Pua. Onecan specifynperm to obtain a list ofnperm randomlypermuted variables. The list is converted to a single ordinary formulaifnperm=1. Ifnperm is omitted, variables are sorted indescending order of the number ofNAs.reformM alsoprints a recommended number of multiple imputations to use, which is aminimum of 5 and the percent of incomplete observations.

Usage

aregImpute(formula, data, subset, n.impute=5, group=NULL,           nk=3, tlinear=TRUE, type=c('pmm','regression','normpmm'),           pmmtype=1, match=c('weighted','closest','kclosest'),           kclosest=3, fweighted=0.2,           curtail=TRUE, constraint=NULL,           boot.method=c('simple', 'approximate bayesian'),           burnin=3, x=FALSE, pr=TRUE, plotTrans=FALSE, tolerance=NULL, B=75)## S3 method for class 'aregImpute'print(x, digits=3, ...)## S3 method for class 'aregImpute'plot(x, nclass=NULL, type=c('ecdf','hist'),     datadensity=c("hist", "none", "rug", "density"),     diagnostics=FALSE, maxn=10, ...)reformM(formula, data, nperm)

Arguments

Details

The sequence of steps used by thearegImpute algorithm is thefollowing.
(1) For each variable containing mNAs where m > 0, initialize theNAs to values from a random sample (without replacement ifa sufficient number of non-missing values exist) of size m from thenon-missing values.
(2) Forburnin+n.impute iterations do the following steps. Thefirstburnin iterations provide a burn-in, and imputations aresaved only from the lastn.impute iterations.
(3) For each variable containing anyNAs, draw a sample withreplacement from the observations in the entire dataset in which thecurrent variable being imputed is non-missing. Fit a flexibleadditive model to predict this target variable while finding theoptimum transformation of it (unless the identitytransformation is forced). Use this fitted flexible model topredict the target variable in all of the original observations.Impute each missing value of the target variable with the observedvalue whose predicted transformed value is closest to the predictedtransformed value of the missing value (ifmatch="closest" andtype="pmm"), or use a draw from a multinomial distribution with probabilities derivedfrom distance weights, ifmatch="weighted" (the default).
(4) After these imputations are computed, use these random drawimputations the next time the curent target variable is used as apredictor of other sometimes-missing variables.

Whenmatch="closest", predictive mean matching does not work wellwhen fewer than 3 variables are used to predict the target variable,because many of the multiple imputations for an observation will beidentical. In the extreme case of one right-hand-side variable andassuming that only monotonic transformations of left and right-sidevariables are allowed, every bootstrap resample will give predictedvalues of the target variable that are monotonically related topredicted values from every other bootstrap resample. The same is truefor Bayesian predicted values. This causes predictive mean matching toalways match on the same donor observation.

When the missingness mechanism for a variable is so systematic that thedistribution of observed values is truncated, predictive mean matchingdoes not work. It will only yield imputed values that are near observedvalues, so intervals in which no values are observed will not bepopulated by imputed values. For this case, the only hope is to makeregression assumptions and use extrapolation. Withtype="regression",aregImpute will use linearextrapolation to obtain a (hopefully) reasonable distribution of imputedvalues. The"regression" option causesaregImpute toimpute missing values by adding a random sample of residuals (withreplacement if there are moreNAs than measured values) on thetransformed scale of the target variable. After random residuals areadded, predicted random draws are obtained on the original untransformedscale using reverse linear interpolation on the table of original andtransformed target values (linear extrapolation when a random residualis large enough to put the random draw prediction outside the range ofobserved values). The bootstrap is used as withtype="pmm" tofactor in the uncertainty of the imputation model.

As model uncertainty is high when the transformation of a targetvariable is unknown,tlinear defaults toTRUE to limit thevariance in predicted values whennk is positive.

Value

a list of class"aregImpute" containing the following elements:

Author(s)

Frank Harrell
Department of Biostatistics
Vanderbilt University
fh@fharrell.com

References

van Buuren, Stef. Flexible Imputation of Missing Data. Chapman &Hall/CRC, Boca Raton FL, 2012.

Little R, An H. Robust likelihood-based analysis of multivariate datawith missing values. Statistica Sinica 14:949-968, 2004.

van Buuren S, Brand JPL, Groothuis-Oudshoorn CGM, Rubin DB. Fullyconditional specifications in multivariate imputation. J Stat CompSim 72:1049-1064, 2006.

de Groot JAH, Janssen KJM, Zwinderman AH, Moons KGM, Reitsma JB.Multiple imputation to correct for partial verification biasrevisited. Stat Med 27:5880-5889, 2008.

Siddique J. Multiple imputation using an iterative hot-deck withdistance-based donor selection. Stat Med 27:83-102, 2008.

White IR, Royston P, Wood AM. Multiple imputation using chainedequations: Issues and guidance for practice. Stat Med 30:377-399,2011.

Curnow E, Carpenter JR, Heron JE, et al: Multiple imputation ofmissing data under missing at random: compatible imputation models arenot sufficient to avoid bias if they are mis-specified. J Clin EpiJune 9, 2023. DOI:10.1016/j.jclinepi.2023.06.011.

See Also

fit.mult.impute,transcan,areg,naclus,naplot,mice,dotchart3,Ecdf,completer

Examples

# Check that aregImpute can almost exactly estimate missing values when# there is a perfect nonlinear relationship between two variables# Fit restricted cubic splines with 4 knots for x1 and x2, linear for x3set.seed(3)x1 <- rnorm(200)x2 <- x1^2x3 <- runif(200)m <- 30x2[1:m] <- NAa <- aregImpute(~x1+x2+I(x3), n.impute=5, nk=4, match='closest')amatplot(x1[1:m]^2, a$imputed$x2)abline(a=0, b=1, lty=2)x1[1:m]^2a$imputed$x2# Multiple imputation and estimation of variances and covariances of# regression coefficient estimates accounting for imputation# Example 1: large sample size, much missing data, no overlap in# NAs across variablesx1 <- factor(sample(c('a','b','c'),1000,TRUE))x2 <- (x1=='b') + 3*(x1=='c') + rnorm(1000,0,2)x3 <- rnorm(1000)y  <- x2 + 1*(x1=='c') + .2*x3 + rnorm(1000,0,2)orig.x1 <- x1[1:250]orig.x2 <- x2[251:350]x1[1:250] <- NAx2[251:350] <- NAd <- data.frame(x1,x2,x3,y, stringsAsFactors=TRUE)# Find value of nk that yields best validating imputation models# tlinear=FALSE means to not force the target variable to be linearf <- aregImpute(~y + x1 + x2 + x3, nk=c(0,3:5), tlinear=FALSE,                data=d, B=10) # normally B=75f# Try forcing target variable (x1, then x2) to be linear while allowing# predictors to be nonlinear (could also say tlinear=TRUE)f <- aregImpute(~y + x1 + x2 + x3, nk=c(0,3:5), data=d, B=10)f## Not run: # Use 100 imputations to better check against individual true valuesf <- aregImpute(~y + x1 + x2 + x3, n.impute=100, data=d)fpar(mfrow=c(2,1))plot(f)modecat <- function(u) { tab <- table(u) as.numeric(names(tab)[tab==max(tab)][1])}table(orig.x1,apply(f$imputed$x1, 1, modecat))par(mfrow=c(1,1))plot(orig.x2, apply(f$imputed$x2, 1, mean))fmi <- fit.mult.impute(y ~ x1 + x2 + x3, lm, f,                        data=d)sqrt(diag(vcov(fmi)))fcc <- lm(y ~ x1 + x2 + x3)summary(fcc)   # SEs are larger than from mult. imputation## End(Not run)## Not run: # Example 2: Very discriminating imputation models,# x1 and x2 have some NAs on the same rows, smaller nset.seed(5)x1 <- factor(sample(c('a','b','c'),100,TRUE))x2 <- (x1=='b') + 3*(x1=='c') + rnorm(100,0,.4)x3 <- rnorm(100)y  <- x2 + 1*(x1=='c') + .2*x3 + rnorm(100,0,.4)orig.x1 <- x1[1:20]orig.x2 <- x2[18:23]x1[1:20] <- NAx2[18:23] <- NA#x2[21:25] <- NAd <- data.frame(x1,x2,x3,y, stringsAsFactors=TRUE)n <- naclus(d)plot(n); naplot(n)  # Show patterns of NAs# 100 imputations to study them; normally use 5 or 10f  <- aregImpute(~y + x1 + x2 + x3, n.impute=100, nk=0, data=d)par(mfrow=c(2,3))plot(f, diagnostics=TRUE, maxn=2)# Note: diagnostics=TRUE makes graphs similar to those made by:# r <- range(f$imputed$x2, orig.x2)# for(i in 1:6) {  # use 1:2 to mimic maxn=2#   plot(1:100, f$imputed$x2[i,], ylim=r,#        ylab=paste("Imputations for Obs.",i))#   abline(h=orig.x2[i],lty=2)# }table(orig.x1,apply(f$imputed$x1, 1, modecat))par(mfrow=c(1,1))plot(orig.x2, apply(f$imputed$x2, 1, mean))fmi <- fit.mult.impute(y ~ x1 + x2, lm, f,                        data=d)sqrt(diag(vcov(fmi)))fcc <- lm(y ~ x1 + x2)summary(fcc)   # SEs are larger than from mult. imputation## End(Not run)## Not run: # Study relationship between smoothing parameter for weighting function# (multiplier of mean absolute distance of transformed predicted# values, used in tricube weighting function) and standard deviation# of multiple imputations.  SDs are computed from average variances# across subjects.  match="closest" same as match="weighted" with# small value of fweighted.# This example also shows problems with predicted mean# matching almost always giving the same imputed values when there is# only one predictor (regression coefficients change over multiple# imputations but predicted values are virtually 1-1 functions of each# other)set.seed(23)x <- runif(200)y <- x + runif(200, -.05, .05)r <- resid(lsfit(x,y))rmse <- sqrt(sum(r^2)/(200-2))   # sqrt of residual MSEy[1:20] <- NAd <- data.frame(x,y)f <- aregImpute(~ x + y, n.impute=10, match='closest', data=d)# As an aside here is how to create a completed dataset for imputation# number 3 as fit.mult.impute would do automatically.  In this degenerate# case changing 3 to 1-2,4-10 will not alter the results.imputed <- impute.transcan(f, imputation=3, data=d, list.out=TRUE,                           pr=FALSE, check=FALSE)sd <- sqrt(mean(apply(f$imputed$y, 1, var)))ss <- c(0, .01, .02, seq(.05, 1, length=20))sds <- ss; sds[1] <- sdfor(i in 2:length(ss)) {  f <- aregImpute(~ x + y, n.impute=10, fweighted=ss[i])  sds[i] <- sqrt(mean(apply(f$imputed$y, 1, var)))}plot(ss, sds, xlab='Smoothing Parameter', ylab='SD of Imputed Values',     type='b')abline(v=.2,  lty=2)  # default value of fweightedabline(h=rmse, lty=2)  # root MSE of residuals from linear regression## End(Not run)## Not run: # Do a similar experiment for the Titanic datasetgetHdata(titanic3)h <- lm(age ~ sex + pclass + survived, data=titanic3)rmse <- summary(h)$sigmaset.seed(21)f <- aregImpute(~ age + sex + pclass + survived, n.impute=10,                data=titanic3, match='closest')sd <- sqrt(mean(apply(f$imputed$age, 1, var)))ss <- c(0, .01, .02, seq(.05, 1, length=20))sds <- ss; sds[1] <- sdfor(i in 2:length(ss)) {  f <- aregImpute(~ age + sex + pclass + survived, data=titanic3,                  n.impute=10, fweighted=ss[i])  sds[i] <- sqrt(mean(apply(f$imputed$age, 1, var)))}plot(ss, sds, xlab='Smoothing Parameter', ylab='SD of Imputed Values',     type='b')abline(v=.2,   lty=2)  # default value of fweightedabline(h=rmse, lty=2)  # root MSE of residuals from linear regression## End(Not run)set.seed(2)d <- data.frame(x1=runif(50), x2=c(rep(NA, 10), runif(40)),                x3=c(runif(4), rep(NA, 11), runif(35)))reformM(~ x1 + x2 + x3, data=d)reformM(~ x1 + x2 + x3, data=d, nperm=2)# Give result or one of the results as the first argument to aregImpute# Constrain imputed values for two variables# Require imputed values for x2 to be above 0.2# Assume x1 is never missing and require imputed values for# x3 to be less than the recipient's value of x1a <- aregImpute(~ x1 + x2 + x3, data=d,                constraint=list(x2 = expression(d$x2 > 0.2),                                x3 = expression(d$x3 < r$x1)))a

Bivariate Summaries Computed Separately by a Series of Predictors

Description

biVar is a generic function that accepts a formula and usualdata,subset, andna.action parameters plus aliststatinfo that specifies a function of two variables tocompute along with information about labeling results for printing andplotting. The function is called separately with each right hand sidevariable and the same left hand variable. The result is a matrix ofbivariate statistics and thestatinfo list that drives printingand plotting. The plot method draws a dot plot with x-axis values bydefault sorted in order of one of the statistics computed by the function.

spearman2 computes the square of Spearman's rho rank correlationand a generalization of it in whichx can relatenon-monotonically toy. This is done by computing the Spearmanmultiple rho-squared between(rank(x), rank(x)^2) andy.Whenx is categorical, a different kind of Spearman correlationused in the Kruskal-Wallis test is computed (andspearman2 can dothe Kruskal-Wallis test). This is done by computing the ordinarymultipleR^2 betweenk-1 dummy variables andrank(y), wherex hask categories.x canalso be a formula, in which case each predictor is correlated separatelywithy, using non-missing observations for that predictor.biVar is used to do the looping and bookkeeping. By default theplot shows the adjustedrho^2, using the same formula used forthe ordinary adjustedR^2. TheF test uses the unadjustedR2.

spearman computes Spearman's rho on non-missing values of twovariables.spearman.test is a simple version ofspearman2.default.

chiSquare is set up likespearman2 except it is intendedfor a categorical response variable. Separate Pearson chi-square testsare done for each predictor, with optional collapsing of infrequentcategories. Numeric predictors having more thang levels arecategorized intog quantile groups.chiSquare usesbiVar.

Usage

biVar(formula, statinfo, data=NULL, subset=NULL,      na.action=na.retain, exclude.imputed=TRUE, ...)## S3 method for class 'biVar'print(x, ...)## S3 method for class 'biVar'plot(x, what=info$defaultwhat,                       sort.=TRUE, main, xlab,                       vnames=c('names','labels'), ...)spearman2(x, ...)## Default S3 method:spearman2(x, y, p=1, minlev=0, na.rm=TRUE, exclude.imputed=na.rm, ...)## S3 method for class 'formula'spearman2(formula, data=NULL,          subset, na.action=na.retain, exclude.imputed=TRUE, ...)spearman(x, y)spearman.test(x, y, p=1)chiSquare(formula, data=NULL, subset=NULL, na.action=na.retain,          exclude.imputed=TRUE, ...)

Arguments

Details

Uses midranks in case of ties, as described by Hollander and Wolfe.P-values for Spearman, Wilcoxon, or Kruskal-Wallis tests areapproximated by using thet orF distributions.

Value

spearman2.default (thefunction that is called for a singlex, i.e., when there is noformula) returns a vector of statistics for the variable.biVar,spearman2.formula, andchiSquare return amatrix with rows corresponding to predictors.

Author(s)

Frank Harrell
Department of Biostatistics
Vanderbilt University
fh@fharrell.com

References

Hollander M. and Wolfe D.A. (1973). Nonparametric Statistical Methods.New York: Wiley.

Press WH, Flannery BP, Teukolsky SA, Vetterling, WT (1988): NumericalRecipes in C. Cambridge: Cambridge University Press.

See Also

combine.levels,varclus,dotchart3,impute,chisq.test,cut2.

Examples

x <- c(-2, -1, 0, 1, 2)y <- c(4,   1, 0, 1, 4)z <- c(1,   2, 3, 4, NA)v <- c(1,   2, 3, 4, 5)spearman2(x, y)plot(spearman2(z ~ x + y + v, p=2))f <- chiSquare(z ~ x + y + v)f

Confidence Intervals for Binomial Probabilities

Description

Produces 1-alpha confidence intervals for binomial probabilities.

Usage

binconf(x, n, alpha=0.05,        method=c("wilson","exact","asymptotic","all"),        include.x=FALSE, include.n=FALSE, return.df=FALSE)

Arguments

Value

a matrix or data.frame containing the computed intervals and,optionally,x andn.

Author(s)

Rollin Brant, Modified by Frank Harrell and
Brad Biggerstaff
Centers for Disease Control and Prevention
National Center for Infectious Diseases
Division of Vector-Borne Infectious Diseases
P.O. Box 2087, Fort Collins, CO, 80522-2087, USA
bkb5@cdc.gov

References

A. Agresti and B.A. Coull, Approximate is better than "exact" forinterval estimation of binomial proportions,American Statistician,52:119–126, 1998.

R.G. Newcombe, Logit confidence intervals and the inverse sinhtransformation,American Statistician,55:200–202, 2001.

L.D. Brown, T.T. Cai and A. DasGupta, Interval estimation fora binomial proportion (with discussion),Statistical Science,16:101–133, 2001.

Examples

binconf(0:10,10,include.x=TRUE,include.n=TRUE)binconf(46,50,method="all")

Bootstrap Kaplan-Meier Estimates

Description

Bootstraps Kaplan-Meier estimate of the probability of survival to atleast a fixed time (times variable) or the estimate of theqquantile of the survival distribution (e.g., median survival time, thedefault).

Usage

bootkm(S, q=0.5, B=500, times, pr=TRUE)

Arguments

Details

bootkm uses Therneau'ssurvfitKM function to efficientlycompute Kaplan-Meier estimates.

Value

a vector containingB bootstrap estimates

Side Effects

updates.Random.seed, and, ifpr=TRUE, prints progressof simulations

Author(s)

Frank Harrell
Department of Biostatistics
Vanderbilt University School of Medicine
fh@fharrell.com

References

Akritas MG (1986): Bootstrapping the Kaplan-Meier estimator. JASA81:1032–1038.

See Also

survfit,Surv,Survival.cph,Quantile.cph

Examples

# Compute 0.95 nonparametric confidence interval for the difference in# median survival time between females and males (two-sample problem)set.seed(1)library(survival)S <- Surv(runif(200))      # no censoringsex <- c(rep('female',100),rep('male',100))med.female <- bootkm(S[sex=='female',], B=100) # normally B=500med.male   <- bootkm(S[sex=='male',],   B=100)describe(med.female-med.male)quantile(med.female-med.male, c(.025,.975), na.rm=TRUE)# na.rm needed because some bootstrap estimates of median survival# time may be missing when a bootstrap sample did not include the# longer survival times

Power and Sample Size for Two-Sample Binomial Test

Description

Uses method of Fleiss, Tytun, and Ury (but without the continuitycorrection) to estimate the power (or the sample size to achieve a givenpower) of a two-sided test for the difference in two proportions. The twosample sizes are allowed to be unequal, but forbsamsize you must specifythe fraction of observations in group 1. For power calculations, oneprobability (p1) must be given, and either the other probability (p2),anodds.ratio, or apercent.reduction must be given. Forbpower orbsamsize, any or all of the arguments may be vectors, in which case theyreturn a vector of powers or sample sizes. All vector arguments must havethe same length.

Givenp1, p2,ballocation uses the method of Brittain and Schlesselmanto compute the optimal fraction of observations to be placed in group 1that either (1) minimize the variance of the difference in two proportions,(2) minimize the variance of the ratio of the two proportions, (3) minimize the variance of the log odds ratio, or(4) maximize the power of the 2-tailed test for differences. For (4)the total sample size must be given, or the fraction optimizingthe power is not returned. The fraction for (3) is one minus the fractionfor (1).

bpower.sim estimates power by simulations, in minimal time. By usingbpower.sim you can see that the formulas without any continuity correctionare quite accurate, and that the power of a continuity-corrected testis significantly lower. That's why no continuity corrections are implementedhere.

Usage

bpower(p1, p2, odds.ratio, percent.reduction,        n, n1, n2, alpha=0.05)bsamsize(p1, p2, fraction=.5, alpha=.05, power=.8)ballocation(p1, p2, n, alpha=.05)bpower.sim(p1, p2, odds.ratio, percent.reduction,            n, n1, n2,            alpha=0.05, nsim=10000)

Arguments

Details

Forbpower.sim, all arguments must be of length one.

Value

forbpower, the power estimate; forbsamsize, a vector containingthe sample sizes in the two groups; forballocation, a vector with4 fractions of observations allocated to group 1, optimizing the fourcriteria mentioned above. Forbpower.sim, a vector with threeelements is returned, corresponding to the simulated power and itslower and upper 0.95 confidence limits.

AUTHOR

Frank Harrell

Department of Biostatistics

Vanderbilt University

fh@fharrell.com

References

Fleiss JL, Tytun A, Ury HK (1980): A simple approximation for calculatingsample sizes for comparing independent proportions. Biometrics 36:343–6.

Brittain E, Schlesselman JJ (1982): Optimal allocation for the comparisonof proportions. Biometrics 38:1003–9.

Gordon I, Watson R (1996): The myth of continuity-corrected sample sizeformulae. Biometrics 52:71–6.

See Also

samplesize.bin,chisq.test,binconf

Examples

bpower(.1, odds.ratio=.9, n=1000, alpha=c(.01,.05))bpower.sim(.1, odds.ratio=.9, n=1000)bsamsize(.1, .05, power=.95)ballocation(.1, .5, n=100)# Plot power vs. n for various odds ratios  (base prob.=.1)n  <- seq(10, 1000, by=10)OR <- seq(.2,.9,by=.1)plot(0, 0, xlim=range(n), ylim=c(0,1), xlab="n", ylab="Power", type="n")for(or in OR) {  lines(n, bpower(.1, odds.ratio=or, n=n))  text(350, bpower(.1, odds.ratio=or, n=350)-.02, format(or))}# Another way to plot the same curves, but letting labcurve do the# work, including labeling each curve at points of maximum separationpow <- lapply(OR, function(or,n)list(x=n,y=bpower(p1=.1,odds.ratio=or,n=n)),              n=n)names(pow) <- format(OR)labcurve(pow, pl=TRUE, xlab='n', ylab='Power')# Contour graph for various probabilities of outcome in the control# group, fixing the odds ratio at .8 ([p2/(1-p2) / p1/(1-p1)] = .8)# n is varied alsop1 <- seq(.01,.99,by=.01)n  <- seq(100,5000,by=250)pow <- outer(p1, n, function(p1,n) bpower(p1, n=n, odds.ratio=.8))# This forms a length(p1)*length(n) matrix of power estimatescontour(p1, n, pow)

Box-percentile plots

Description

Producess side-by-side box-percentile plots from several vectors or alist of vectors.

Usage

bpplot(..., name=TRUE, main="Box-Percentile Plot",        xlab="", ylab="", srtx=0, plotopts=NULL)

Arguments

Value

There are no returned values

Side Effects

A plot is created on the current graphics device.

BACKGROUND

Box-percentile plots are similiar to boxplots, except box-percentile plotssupply more information about the univariate distributions. At any heightthe width of the irregular "box" is proportional to the percentile of thatheight, up to the 50th percentile, and above the 50th percentile the widthis proportional to 100 minus the percentile. Thus, the width at any givenheight is proportional to the percent of observations that are more extreme in that direction. As in boxplots, the median, 25th and 75th percentiles are marked with line segments across the box.

Author(s)

Jeffrey Banfield
umsfjban@bill.oscs.montana.edu
Modified by F. Harrell 30Jun97

References

Esty WW, Banfield J: The box-percentile plot. J StatisticalSoftware 8 No. 17, 2003.

See Also

panel.bpplot,boxplot,Ecdf,bwplot

Examples

set.seed(1)x1 <- rnorm(500)x2 <- runif(500, -2, 2)x3 <- abs(rnorm(500))-2bpplot(x1, x2, x3)g <- sample(1:2, 500, replace=TRUE)bpplot(split(x2, g), name=c('Group 1','Group 2'))rm(x1,x2,x3,g)

Statistics by Categories

Description

For any number of cross-classification variables,bystatsreturns a matrix with the sample size, number missingy, andfun(non-missing y), with the cross-classifications designatedby rows. Uses Harrell's modification of theinteractionfunction to produce cross-classifications. The defaultfun ismean, and ify is binary, the mean is labeled asFraction. There is aprint method as well as alatex method for objects created bybystats.bystats2 handles the special case in which there are 2classifcation variables, and places the first one in rows and thesecond in columns. Theprint method forbystats2 usestheprint.char.matrix function to organize statisticsfor cells into boxes.

Usage

bystats(y, ..., fun, nmiss, subset)## S3 method for class 'bystats'print(x, ...)## S3 method for class 'bystats'latex(object, title, caption, rowlabel, ...)bystats2(y, v, h, fun, nmiss, subset)## S3 method for class 'bystats2'print(x, abbreviate.dimnames=FALSE,   prefix.width=max(nchar(dimnames(x)[[1]])), ...)## S3 method for class 'bystats2'latex(object, title, caption, rowlabel, ...)

Arguments

Value

forbystats, a matrix with row names equal to the classification labels and columnnamesN, Missing, funlab, wherefunlab is determined fromfun.A row is added to the end with the summary statistics computed on all observations combined. The class of this matrix isbystats.Forbystats, returns a 3-dimensional array with the last dimensioncorresponding to statistics being computed. The class of the array isbystats2.

Side Effects

latex produces a.tex file.

Author(s)

Frank Harrell
Department of Biostatistics
Vanderbilt University
fh@fharrell.com

See Also

interaction,cut,cut2,latex,print.char.matrix,translate

Examples

## Not run: bystats(sex==2, county, city)bystats(death, race)bystats(death, cut2(age,g=5), race)bystats(cholesterol, cut2(age,g=4), sex, fun=median)bystats(cholesterol, sex, fun=quantile)bystats(cholesterol, sex, fun=function(x)c(Mean=mean(x),Median=median(x)))latex(bystats(death,race,nmiss=FALSE,subset=sex=="female"), digits=2)f <- function(y) c(Hazard=sum(y[,2])/sum(y[,1]))# f() gets the hazard estimate for right-censored data from exponential dist.bystats(cbind(d.time, death), race, sex, fun=f)bystats(cbind(pressure, cholesterol), age.decile,         fun=function(y) c(Median.pressure   =median(y[,1]),                          Median.cholesterol=median(y[,2])))y <- cbind(pressure, cholesterol)bystats(y, age.decile,         fun=function(y) apply(y, 2, median))   # same result as last onebystats(y, age.decile, fun=function(y) apply(y, 2, quantile, c(.25,.75)))# The last one computes separately the 0.25 and 0.75 quantiles of 2 vars.latex(bystats2(death, race, sex, fun=table))## End(Not run)

capitalize the first letter of a string

Description

Capitalizes the first letter of each element of the string vector.

Usage

capitalize(string)

Arguments

Value

Returns a vector of charaters with the first letter capitalized

Author(s)

Charles Dupont

Examples

capitalize(c("Hello", "bob", "daN"))

Power of Interaction Test for Exponential Survival

Description

Uses the method of Peterson and George to compute the power of aninteraction test in a 2 x 2 setup in which all 4 distributions areexponential. This will be the same as the power of the Cox modeltest if assumptions hold. The test is 2-tailed. The duration of accrual is specified(constant accrual is assumed), as is the minimum follow-up time.The maximum follow-up time is thenaccrual + tmin. Treatmentallocation is assumed to be 1:1.

Usage

ciapower(tref, n1, n2, m1c, m2c, r1, r2, accrual, tmin,          alpha=0.05, pr=TRUE)

Arguments

Value

power

Side Effects

prints

AUTHOR

Frank Harrell

Department of Biostatistics

Vanderbilt University

fh@fharrell.com

References

Peterson B, George SL: Controlled Clinical Trials 14:511–522; 1993.

See Also

cpower,spower

Examples

# Find the power of a race x treatment test.  25% of patients will# be non-white and the total sample size is 14000.  # Accrual is for 1.5 years and minimum follow-up is 5y.# Reduction in 5-year mortality is 15% for whites, 0% or -5% for# non-whites.  5-year mortality for control subjects if assumed to# be 0.18 for whites, 0.23 for non-whites.n <- 14000for(nonwhite.reduction in c(0,-5)) {  cat("\n\n\n% Reduction in 5-year mortality for non-whites:",      nonwhite.reduction, "\n\n")  pow <- ciapower(5,  .75*n, .25*n,  .18, .23,  15, nonwhite.reduction,                    1.5, 5)  cat("\n\nPower:",format(pow),"\n")}

Convert between the 5 different coordinate sytems on a graphical device

Description

Takes a set of coordinates in any of the 5 coordinate systems (usr,plt, fig, dev, or tdev) and returns the same points in all 5coordinate systems.

Usage

cnvrt.coords(x, y = NULL, input = c("usr", "plt", "fig", "dev","tdev"))

Arguments

Details

Every plot has 5 coordinate systems:

usr (User): the coordinate system of the data, this is shown by thetick marks and axis labels.

plt (Plot): Plot area, coordinates range from 0 to 1 with 0corresponding to the x and y axes and 1 corresponding to the top andright of the plot area. Margins of the plot correspond to plotcoordinates less than 0 or greater than 1.

fig (Figure): Figure area, coordinates range from 0 to 1 with 0corresponding to the bottom and left edges of the figure (includingmargins, label areas) and 1 corresponds to the top and right edges.fig and dev coordinates will be identical if there is only 1 figurearea on the device (layout, mfrow, or mfcol has not been used).

dev (Device): Device area, coordinates range from 0 to 1 with 0corresponding to the bottom and left of the device region within theouter margins and 1 is the top and right of the region withing theouter margins. If the outer margins are all set to 0 then tdev anddev should be identical.

tdev (Total Device): Total Device area, coordinates range from 0 to 1 with 0corresponding to the bottom and left edges of the device (piece ofpaper, window on screen) and 1 corresponds to the top and right edges.

Value

A list with 5 components, each component is a list with vectors namedx and y. The 5 sublists are:

Note

You must provide both x and y, but one of them may beNA.

This function is becoming depricated with the new functionsgrconvertX andgrconvertY in R version 2.7.0 and beyond.These new functions use the correct coordinate system names and havemore coordinate systems available, you should start using them instead.

Author(s)

Greg Snowgreg.snow@imail.org

See Also

par specifically 'usr','plt', and 'fig'. Also'xpd' for plotting outside of the plotting region and 'mfrow' and'mfcol' for multi figure plotting.subplot,grconvertX andgrconvertY in R2.7.0 and later

Examples

old.par <- par(no.readonly=TRUE)par(mfrow=c(2,2),xpd=NA)# generate some sample datatmp.x <- rnorm(25, 10, 2)tmp.y <- rnorm(25, 50, 10)tmp.z <- rnorm(25, 0, 1)plot( tmp.x, tmp.y)# draw a diagonal line across the plot areatmp1 <- cnvrt.coords( c(0,1), c(0,1), input='plt' )lines(tmp1$usr, col='blue')# draw a diagonal line accross figure regiontmp2 <- cnvrt.coords( c(0,1), c(1,0), input='fig')lines(tmp2$usr, col='red')# save coordinate of point 1 and y value near top of plot for future plotstmp.point1 <- cnvrt.coords(tmp.x[1], tmp.y[1])tmp.range1 <- cnvrt.coords(NA, 0.98, input='plt')# make a second plot and draw a line linking point 1 in each plotplot(tmp.y, tmp.z)tmp.point2 <- cnvrt.coords( tmp.point1$dev, input='dev' )arrows( tmp.y[1], tmp.z[1], tmp.point2$usr$x, tmp.point2$usr$y, col='green')# draw another plot and add rectangle showing same range in 2 plotsplot(tmp.x, tmp.z)tmp.range2 <- cnvrt.coords(NA, 0.02, input='plt')tmp.range3 <- cnvrt.coords(NA, tmp.range1$dev$y, input='dev')rect( 9, tmp.range2$usr$y, 11, tmp.range3$usr$y, border='yellow')# put a label just to the right of the plot and#  near the top of the figure region.text( cnvrt.coords(1.05, NA, input='plt')$usr$x,cnvrt.coords(NA, 0.75, input='fig')$usr$y,"Label", adj=0)par(mfrow=c(1,1))## create a subplot within another plot (see also subplot)plot(1:10, 1:10)tmp <- cnvrt.coords( c( 1, 4, 6, 9), c(6, 9, 1, 4) )par(plt = c(tmp$dev$x[1:2], tmp$dev$y[1:2]), new=TRUE)hist(rnorm(100))par(fig = c(tmp$dev$x[3:4], tmp$dev$y[3:4]), new=TRUE)hist(rnorm(100))par(old.par)

Miscellaneous ggplot2 and grid Helper Functions

Description

These functions are used onggplot2 objects or as layers whenbuilding aggplot2 object, and to facilitate use ofgridExtra.colorFacet colors the thin rectangles used to separate panels created byfacet_grid (andprobably byfacet_wrap). A better approach may be found athttps://stackoverflow.com/questions/28652284/.arrGrob is a front-end togridExtra::arrangeGrob thatallows for proper printing. Seehttps://stackoverflow.com/questions/29062766/store-output-from-gridextragrid-arrange-into-an-object/. ThearrGrobprint method callsgrid::grid.draw.

Usage

colorFacet(g, col = adjustcolor("blue", alpha.f = 0.3))arrGrob(...)## S3 method for class 'arrGrob'print(x, ...)

Arguments

Author(s)

Sandy Muspratt

Examples

## Not run: require(ggplot2)s <- summaryP(age + sex ~ region + treatment)colorFacet(ggplot(s))   # prints directly# arrGrob is called by rms::ggplot.Predict and others## End(Not run)

combine.levels

Description

Combine Infrequent Levels of a Categorical Variable

Usage

combine.levels(  x,  minlev = 0.05,  m,  ord = is.ordered(x),  plevels = FALSE,  sep = ",")

Arguments

Details

After turning 'x' into a 'factor' if it is not one already, combineslevels of 'x' whose frequency falls below a specified relative frequency 'minlev' or absolute count 'm'. When 'x' is not treated as ordered, all of thesmall frequency levels are combined into '"OTHER"', unless 'plevels=TRUE'.When 'ord=TRUE' or 'x' is an ordered factor, only consecutive levelsare combined. New levels are constructed by concatenating the levels with'sep' as a separator. This is useful when comparing ordinal regressionwith polytomous (multinomial) regression and there are too manycategories for polytomous regression. 'combine.levels' is also usefulwhen assumptions of ordinal models are being checked empirically bycomputing exceedance probabilities for various cutoffs of thedependent variable.

Value

a factor variable, or if 'ord=TRUE' an ordered factor variable

Author(s)

Frank Harrell

Examples

x <- c(rep('A', 1), rep('B', 3), rep('C', 4), rep('D',1), rep('E',1))combine.levels(x, m=3)combine.levels(x, m=3, plevels=TRUE)combine.levels(x, ord=TRUE, m=3)x <- c(rep('A', 1), rep('B', 3), rep('C', 4), rep('D',1), rep('E',1),       rep('F',1))combine.levels(x, ord=TRUE, m=3)

Combination Plot

Description

Generates a plotly attribute plot given a series of possibly overlapping binary variables

Usage

combplotp(  formula,  data = NULL,  subset,  na.action = na.retain,  vnames = c("labels", "names"),  includenone = FALSE,  showno = FALSE,  maxcomb = NULL,  minfreq = NULL,  N = NULL,  pos = function(x) 1 * (tolower(x) %in% c("true", "yes", "y", "positive", "+",    "present", "1")),  obsname = "subjects",  ptsize = 35,  width = NULL,  height = NULL,  ...)

Arguments

Details

Similar to theUpSetR package, draws a special dot chart sometimes called an attribute plot that depicts all possible combination of the binary variables. By default a positive value, indicating that a certain condition pertains for a subject, is any of logicalTRUE, numeric 1,"yes","y","positive","+" or"present" value, and all others are considered negative. The user can override this determination by specifying her ownpos function. Case is ignored in the variable values.

The plot uses solid dots arranged in a vertical line to indicate which combination of conditions is being considered. Frequencies of all possible combinations are shown above the dot chart. Marginal frequencies of positive values for the input variables are shown to the left of the dot chart. More information for all three of these component symbols is provided in hover text.

Variables are sorted in descending order of marginal frqeuencies and likewise for combinations of variables.

Value

plotly object

Author(s)

Frank Harrell

Examples

if (requireNamespace("plotly")) {  g <- function() sample(0:1, n, prob=c(1 - p, p), replace=TRUE)  set.seed(2); n <- 100; p <- 0.5  x1 <- g(); label(x1) <- 'A long label for x1 that describes it'  x2 <- g()  x3 <- g(); label(x3) <- 'This is<br>a label for x3'  x4 <- g()  combplotp(~ x1 + x2 + x3 + x4, showno=TRUE, includenone=TRUE)  n <- 1500; p <- 0.05  pain       <- g()  anxiety    <- g()  depression <- g()  soreness   <- g()  numbness   <- g()  tiredness  <- g()  sleepiness <- g()  combplotp(~ pain + anxiety + depression + soreness + numbness +            tiredness + sleepiness, showno=TRUE)}

completer

Description

Create imputed dataset(s) usingtranscan andaregImpute objects

Usage

completer(a, nimpute, oneimpute = FALSE, mydata)

Arguments

Details

Similar in function tomice::complete, this function usestranscan andaregImpute objects to impute missing dataand returns the completed dataset(s) as a dataframe or a list.It assumes thattranscan is used for single regression imputation.

Value

      A single or a list of completed dataset(s).

Author(s)

      Yong-Hao Pua, Singapore General Hospital

Examples

## Not run: mtcars$hp[1:5]    <- NAmtcars$wt[1:10]   <- NAmyrform <- ~ wt + hp + I(carb)mytranscan  <- transcan( myrform,  data = mtcars, imputed = TRUE,  pl = FALSE, pr = FALSE, trantab = TRUE, long = TRUE)myareg      <- aregImpute(myrform, data = mtcars, x=TRUE, n.impute = 5)completer(mytranscan)                    # single completed datasetcompleter(myareg, 3, oneimpute = TRUE)# single completed dataset based on the `n.impute`th set of multiple imputationcompleter(myareg, 3)# list of completed datasets based on first `nimpute` sets of multiple imputationcompleter(myareg)# list of completed datasets based on all available sets of multiple imputation# To get a stacked data frame of all completed datasets use# do.call(rbind, completer(myareg, data=mydata))# or use rbindlist in data.table## End(Not run)

Element Merging

Description

Merges an object by the names of its elements. Inserting elements invalue intox that do not exists inx andreplacing elements inx that exists invalue withvalue elements ifprotect is false.

Usage

consolidate(x, value, protect, ...)## Default S3 method:consolidate(x, value, protect=FALSE, ...)consolidate(x, protect, ...) <- value

Arguments

Author(s)

Charles Dupont

See Also

names

Examples

x <- 1:5names(x) <- LETTERS[x]y <- 6:10names(y) <- LETTERS[y-2]x                  # c(A=1,B=2,C=3,D=4,E=5)y                  # c(D=6,E=7,F=8,G=9,H=10)consolidate(x, y)      # c(A=1,B=2,C=3,D=6,E=7,F=8,G=9,H=10)consolidate(x, y, protect=TRUE)      # c(A=1,B=2,C=3,D=4,E=5,F=8,G=9,H=10)

Metadata for a Data Frame

Description

contents is a generic method for whichcontents.data.frameis currently the only method.contents.data.frame creates anobject containing the following attributes of the variables from a data frame: names, labels (if any), units (if any), number offactor levels (if any), factor levels,class, storage mode, and number of NAs.print.contents.data.framewill print the results, with options for sorting the variables.html.contents.data.frame creates HTML code for displaying theresults. This code has hyperlinks so that if the user clicks on thenumber of levels the browser jumps to the correct part of a table offactor levels for all thefactor variables. If long labels arepresent ("longlabel" attributes on variables), these are printedat the bottom and thehtml method links to them through theregular labels. Variables having the samelevels in the sameorder have the levels factored out for brevity.

contents.list prints a directory of datasets whensasxport.get imported more than one SAS dataset.

Ifoptions(prType='html') is in effect, callingprint onan object that is the contents of a data frame will result inrendering the HTML version. If run from the console a browser windowwill open.

Usage

contents(object, ...)## S3 method for class 'data.frame'contents(object, sortlevels=FALSE, id=NULL,  range=NULL, values=NULL, ...)## S3 method for class 'contents.data.frame'print(x,    sort=c('none','names','labels','NAs'), prlevels=TRUE, maxlevels=Inf,    number=FALSE, ...) ## S3 method for class 'contents.data.frame'html(object,           sort=c('none','names','labels','NAs'), prlevels=TRUE, maxlevels=Inf,           levelType=c('list','table'),           number=FALSE, nshow=TRUE, ...)## S3 method for class 'list'contents(object, dslabels, ...)## S3 method for class 'contents.list'print(x,    sort=c('none','names','labels','NAs','vars'), ...)

Arguments

Value

an object of class"contents.data.frame" or"contents.list". For thehtml method is anhtmlcharacter vector object.

Author(s)

Frank Harrell
Vanderbilt University
fh@fharrell.com

See Also

describe,html,upData,extractlabs,hlab

Examples

set.seed(1)dfr <- data.frame(x=rnorm(400),y=sample(c('male','female'),400,TRUE),                  stringsAsFactors=TRUE)contents(dfr)dfr <- upData(dfr, labels=c(x='Label for x', y='Label for y'))attr(dfr$x, 'longlabel') <- 'A very long label for x that can continue onto multiple long lines of text'k <- contents(dfr)print(k, sort='names', prlevels=FALSE)## Not run: html(k)html(contents(dfr))            # same resultlatex(k$contents)              # latex.default just the main information## End(Not run)

Power of Cox/log-rank Two-Sample Test

Description

Assumes exponential distributions for both treatment groups.Uses the George-Desu method along withformulas of Schoenfeld that allow estimation of the expected number ofevents in the two groups. To allow for drop-ins (noncompliance to control therapy, crossover tointervention) and noncompliance of the intervention, the method ofLachin and Foulkes is used.

Usage

cpower(tref, n, mc, r, accrual, tmin, noncomp.c=0, noncomp.i=0,        alpha=0.05, nc, ni, pr=TRUE)

Arguments

Details

For handling noncompliance, uses a modification of formula (5.4) ofLachin and Foulkes. Their method is based on a test for the differencein two hazard rates, whereascpower is based on testing the differencein two log hazards. It is assumed here that the same correction factorcan be approximately applied to the log hazard ratio as Lachin and Foulkes applied tothe hazard difference.

Note that Schoenfeld approximates the varianceof the log hazard ratio by4/m, wherem is the total number of events,whereas the George-Desu method uses the slightly better1/m1 + 1/m2.Power from this function will thus differ slightly from that obtained withthe SASsamsizc program.

Value

power

Side Effects

prints

Author(s)

Frank Harrell
Department of Biostatistics
Vanderbilt University
fh@fharrell.com

References

Peterson B, George SL: Controlled Clinical Trials 14:511–522; 1993.

Lachin JM, Foulkes MA: Biometrics 42:507–519; 1986.

Schoenfeld D: Biometrics 39:499–503; 1983.

See Also

spower,ciapower,bpower

Examples

#In this example, 4 plots are drawn on one page, one plot for each#combination of noncompliance percentage.  Within a plot, the#5-year mortality % in the control group is on the x-axis, and#separate curves are drawn for several % reductions in mortality#with the intervention.  The accrual period is 1.5y, with all#patients followed at least 5y and some 6.5y.par(mfrow=c(2,2),oma=c(3,0,3,0))morts <- seq(10,25,length=50)red <- c(10,15,20,25)for(noncomp in c(0,10,15,-1)) {  if(noncomp>=0) nc.i <- nc.c <- noncomp else {nc.i <- 25; nc.c <- 15}  z <- paste("Drop-in ",nc.c,"%, Non-adherence ",nc.i,"%",sep="")  plot(0,0,xlim=range(morts),ylim=c(0,1),           xlab="5-year Mortality in Control Patients (%)",           ylab="Power",type="n")  title(z)  cat(z,"\n")  lty <- 0  for(r in red) {        lty <- lty+1        power <- morts        i <- 0        for(m in morts) {          i <- i+1          power[i] <- cpower(5, 14000, m/100, r, 1.5, 5, nc.c, nc.i, pr=FALSE)        }        lines(morts, power, lty=lty)  }  if(noncomp==0)legend(18,.55,rev(paste(red,"% reduction",sep="")),           lty=4:1,bty="n")}mtitle("Power vs Non-Adherence for Main Comparison",           ll="alpha=.05, 2-tailed, Total N=14000",cex.l=.8)## Point sample size requirement vs. mortality reduction# Root finder (uniroot()) assumes needed sample size is between# 1000 and 40000#nc.i <- 25; nc.c <- 15; mort <- .18red <- seq(10,25,by=.25)samsiz <- redi <- 0for(r in red) {  i <- i+1  samsiz[i] <- uniroot(function(x) cpower(5, x, mort, r, 1.5, 5,                                          nc.c, nc.i, pr=FALSE) - .8,                       c(1000,40000))$root}samsiz <- samsiz/1000par(mfrow=c(1,1))plot(red, samsiz, xlab='% Reduction in 5-Year Mortality', ylab='Total Sample Size (Thousands)', type='n')lines(red, samsiz, lwd=2)title('Sample Size for Power=0.80\nDrop-in 15%, Non-adherence 25%')title(sub='alpha=0.05, 2-tailed', adj=0)

Read Comma-Separated Text Data Files

Description

Read comma-separated text data files, allowing optional translationto lower case for variable names after making them valid S names.There is a facility for reading long variable labels as one of therows. If labels are not specified and a final variable name is notthe same as that in the header, the original variable name is saved asa variable label. Usesread.csv if thedata.tablepackage is not in effect, otherwise callsfread.

Usage

csv.get(file, lowernames=FALSE, datevars=NULL, datetimevars=NULL,        dateformat='%F',        fixdates=c('none','year'), comment.char="", autodate=TRUE,        allow=NULL, charfactor=FALSE,        sep=',', skip=0, vnames=NULL, labels=NULL, text=NULL, ...)

Arguments

Details

csv.get reads comma-separated text data files, allowing optionaltranslation to lower case for variable names after making them valid Snames. Original possibly non-legal names are taken to be variablelabels iflabels is not specified. Character or factorvariables containing dates can be converted to date variables.cleanup.import is invoked to finish the job.

Value

a new data frame.

Author(s)

Frank Harrell, Vanderbilt University

See Also

sas.get,data.frame,cleanup.import,read.csv,strptime,POSIXct,Date,fread

Examples

## Not run: dat <- csv.get('myfile.csv')# Read a csv file with junk in the first row, variable names in the# second, long variable labels in the third, and junk in the 4th rowdat <- csv.get('myfile.csv', vnames=2, labels=3, skip=4)## End(Not run)

Representative Curves

Description

curveRep finds representative curves from arelatively large collection of curves. The curves usually representtime-response profiles as in serial (longitudinal or repeated) datawith possibly unequal time points and greatly varying sample sizes persubject. After excluding records containing missingx ory, records are first stratified intokn groups having similarsample sizes per curve (subject). Within these strata, curves arenext stratified according to the distribution ofx points percurve (typically measurement times per subject). Theclara clustering/partitioning function is usedto do this, clustering on one, two, or threex characteristicsdepending on the minimum sample size in the current interval of samplesize. If the interval has a minimum number of uniquevalues ofone, clustering is done on the singlex values. If the minimumnumber of uniquex values is two, clustering is done to creategroups that are similar on bothmin(x) andmax(x). Forgroups containing no fewer than three uniquex values,clustering is done on the trio of valuesmin(x),max(x),and the longest gap between any successivex. Then withinsample size andx distribution strata, clustering oftime-response profiles is based onp values ofy allevaluated at the samep equally-spacedx's within thestratum. An option allows per-curve data to be smoothed withlowess before proceeding. Outerx values aretaken as extremes ofx across all curves within the stratum.Linear interpolation within curves is used to estimatey at thegrid ofx's. For curves within the stratum that do not extendto the most extremex values in that stratum, extrapolationuses flat lines from the observed extremes in the curve unlessextrap=TRUE. Thepy values are clustered usingclara.

print andplot methods show results. By specifying anauxiliaryidcol variable toplot, other variables suchas treatment may be depicted to allow the analyst to determine forexample whether subjects on different treatments are assigned todifferent time-response profiles. To write the frequencies of avariable such as treatment in the upper left corner of each panel(instead of the grand total number of clusters in that panel), specifyfreq.

curveSmooth takes a set of curves and smooths them usinglowess. If the number of uniquex points in a curve isless thanp, the smooth is evaluated at the uniquexvalues. Otherwise it is evaluated at an equally spaced set ofx points over the observed range. If fewer than 3 uniquex values are in a curve, those points are used and smoothing is not done.

Usage

curveRep(x, y, id, kn = 5, kxdist = 5, k = 5, p = 5,         force1 = TRUE, metric = c("euclidean", "manhattan"),         smooth=FALSE, extrap=FALSE, pr=FALSE)## S3 method for class 'curveRep'print(x, ...)## S3 method for class 'curveRep'plot(x, which=1:length(res),                        method=c('all','lattice','data'),                        m=NULL, probs=c(.5, .25, .75), nx=NULL, fill=TRUE,                        idcol=NULL, freq=NULL, plotfreq=FALSE,                        xlim=range(x), ylim=range(y),                        xlab='x', ylab='y', colorfreq=FALSE, ...)curveSmooth(x, y, id, p=NULL, pr=TRUE)

Arguments

Details

In the graph titles for the default graphic output,n refers to theminimum sample size,x refers to the sequential x-distributioncluster, andc refers to the sequential x-y profile cluster. Graphsfrommethod = "lattice" are produced byxyplot and in the panel titlesdistribution refers to the x-distribution stratum andcluster refers to the x-y profile cluster.

Value

a list of class"curveRep" with the following elements

curveSmooth returns a list with elementsx,y,id.

Note

The references describe other methods for derivingrepresentative curves, but those methods were not used here. The lastreference which used a cluster analysis on principal componentsmotivatedcurveRep however. Thekml package does k-means clustering of longitudinal data with imputation.

Author(s)

Frank Harrell
Department of Biostatistics
Vanderbilt University
fh@fharrell.com

References

Segal M. (1994): Representative curves for longitudinal data viaregression trees. J Comp Graph Stat 3:214-233.

Jones MC, Rice JA (1992): Displaying the important features of largecollections of similar curves. Am Statistician 46:140-145.

Zheng X, Simpson JA, et al (2005): Data from a study of effectivenesssuggested potential prognostic factors related to the patterns ofshoulder pain. J Clin Epi 58:823-830.

See Also

clara,dataRep

Examples

## Not run: # Simulate 200 curves with per-curve sample sizes ranging from 1 to 10# Make curves with odd-numbered IDs have an x-distribution that is random# uniform [0,1] and those with even-numbered IDs have an x-dist. that is# half as wide but still centered at 0.5.  Shift y values higher with# increasing IDsset.seed(1)N <- 200nc <- sample(1:10, N, TRUE)id <- rep(1:N, nc)x <- y <- idfor(i in 1:N) {  x[id==i] <- if(i %% 2) runif(nc[i]) else runif(nc[i], c(.25, .75))  y[id==i] <- i + 10*(x[id==i] - .5) + runif(nc[i], -10, 10)}w <- curveRep(x, y, id, kxdist=2, p=10)wpar(ask=TRUE, mfrow=c(4,5))plot(w)                # show everything, profiles going acrosspar(mfrow=c(2,5))plot(w,1)              # show n=1 results# Use a color assignment table, assigning low curves to green and# high to red.  Unique curve (subject) IDs are the names of the vector.cols <- c(rep('green', N/2), rep('red', N/2))names(cols) <- as.character(1:N)plot(w, 3, idcol=cols)par(ask=FALSE, mfrow=c(1,1))plot(w, 1, 'lattice')  # show n=1 resultsplot(w, 3, 'lattice')  # show n=4-5 resultsplot(w, 3, 'lattice', idcol=cols)  # same but different color mappingplot(w, 3, 'lattice', m=1)  # show a single "representative" curve# Show median, 10th, and 90th percentiles of supposedly representative curvesplot(w, 3, 'lattice', m='quantiles', probs=c(.5,.1,.9))# Same plot but with much less grouping of x variableplot(w, 3, 'lattice', m='quantiles', probs=c(.5,.1,.9), nx=2)# Use ggplot2 for one sample size intervalz <- plot(w, 2, 'data')require(ggplot2)ggplot(z, aes(x, y, color=curve)) + geom_line() +       facet_grid(distribution ~ cluster) +       theme(legend.position='none') +       labs(caption=z$ninterval[1])# Smooth data before profiling.  This allows later plotting to plot# smoothed representative curves rather than raw curves (which# specifying smooth=TRUE to curveRep would do, if curveSmooth was not used)d <- curveSmooth(x, y, id)w <- with(d, curveRep(x, y, id))# Example to show that curveRep can cluster profiles correctly when# there is no noise.  In the data there are four profiles - flat, flat# at a higher mean y, linearly increasing then flat, and flat at the# first height except for a sharp triangular peakset.seed(1)x <- 0:100m <- length(x)profile <- matrix(NA, nrow=m, ncol=4)profile[,1] <- rep(0, m)profile[,2] <- rep(3, m)profile[,3] <- c(0:3, rep(3, m-4))profile[,4] <- c(0,1,3,1,rep(0,m-4))col <- c('black','blue','green','red')matplot(x, profile, type='l', col=col)xeval <- seq(0, 100, length.out=5)s <- x matplot(x[s], profile[s,], type='l', col=col)id <- rep(1:100, each=m)X <- Y <- idcols <- character(100)names(cols) <- as.character(1:100)for(i in 1:100) {  s <- id==i  X[s] <- x  j <- sample(1:4,1)  Y[s] <- profile[,j]  cols[i] <- col[j]}table(cols)yl <- c(-1,4)w <- curveRep(X, Y, id, kn=1, kxdist=1, k=4)plot(w, 1, 'lattice', idcol=cols, ylim=yl)# Found 4 clusters but two have same profilew <- curveRep(X, Y, id, kn=1, kxdist=1, k=3)plot(w, 1, 'lattice', idcol=cols, freq=cols, plotfreq=TRUE, ylim=yl)# Incorrectly combined black and red because default value p=5 did# not result in different profiles at x=xevalw <- curveRep(X, Y, id, kn=1, kxdist=1, k=4, p=40)plot(w, 1, 'lattice', idcol=cols, ylim=yl)# Found correct clusters because evaluated curves at 40 equally# spaced points and could find the sharp triangular peak in profile 4## End(Not run)

Cut a Numeric Variable into Intervals

Description

cut2 is a function likecut but left endpoints are inclusive and labels are ofthe form[lower, upper), except that last interval is[lower,upper]. If cuts are given, will by default make sure that cuts include entirerange ofx.Also, if cuts are not given, will cutx into quantile groups (g given) or groupswith a given minimum number of observations (m). Whereas cut creates acategory object,cut2 creates a factor object.m is not guaranteed but is a target.

cutGn guarantees that the grouped variable will have a minimum ofm observations in any group. This is done by an exhaustive algorithm that runs fast due to being coded in Fortran.

Usage

cut2(x, cuts, m=150, g, levels.mean=FALSE, digits, minmax=TRUE,oneval=TRUE, onlycuts=FALSE, formatfun=format, ...)cutGn(x, m, what=c('mean', 'factor', 'summary', 'cuts', 'function'), rcode=FALSE)

Arguments

Value

a factor variable with levels of the form[a,b) or formatted means(character strings) unlessonlycuts isTRUE in which casea numeric vector is returned

See Also

cut,quantile,combine.levels

Examples

set.seed(1)x <- runif(1000, 0, 100)z <- cut2(x, c(10,20,30))table(z)table(cut2(x, g=10))      # quantile groupstable(cut2(x, m=50))      # group x into intevals with at least 50 obs.table(cutGn(x, m=50, what='factor'))f <- cutGn(x, m=50, what='function')ff(c(-1, 2, 10), what='mean')f(c(-1, 2, 10), what='factor')## Not run:   x <- round(runif(200000), 3)  system.time(a <- cutGn(x, m=20))              # 0.02s  system.time(b <- cutGn(x, m=20, rcode=TRUE))  # 1.51s  identical(a, b)## End(Not run)

Tips for Creating, Modifying, and Checking Data Frames

Description

This help file contains a template for importing data to create an Rdata frame, correcting some problems resulting from the import andmaking the data frame be stored more efficiently, modifying the dataframe (including better annotating it and changing the names of someof its variables), and checking and inspecting the data frame forreasonableness of the values of its variables and to describe patternsof missing data. Various built-in functions and functions in theHmisc library are used. At the end some methods for creating dataframes “from scratch” withinR are presented.

The examples below attempt to clarify the separation of operationsthat are done on a data frame as a whole, operations that are done ona small subset of its variables without attaching the whole dataframe, and operations that are done on many variables after attachingthe data frame in search position one. It also tries to clarify thatfor analyzing several separate variables usingR commands that do notsupport adata argument, it is helpful to attach the data framein a search position later than position one.

It is often useful to create, modify, and process datasets in thefollowing order.

1. Import external data into a data frame (if the raw data do notcontain column names, provide these during the import if possible)

2. Make global changes to a data frame (e.g., changing variablenames)

3. Change attributes or values of variables within a data frame

4. Do analyses involving the whole data frame (without attaching it)
   (Data frame still in .Data)

5. Do analyses of individual variables (after attaching the dataframe in search position two or later)

Details

The examples below use theFEV dataset fromRosner 1995. Almost any dataset would do. The jcetable dataare taken fromGalobardes, etal.

Presently, giving a variable the"units" attribute (using theHmiscunits function) only benefits theHmiscdescribe function and thermslibrary's version of thelink[rms]{Surv} function. Variableslabels defined with the Hmisclabel function are used bydescribe,summary.formula, and many ofthe plotting functions inHmisc andrms.

References

Alzola CF, Harrell FE (2006):An Introduction to S and the Hmisc and Design Libraries.Chapters 3 and 4,https://hbiostat.org/R/doc/sintro.pdf.

Galobardes, et al. (1998),J Clin Epi 51:875-881.

Rosner B (1995):Fundamentals of Biostatistics, 4th Edition.New York: Duxbury Press.

See Also

scan,read.table,cleanup.import,sas.get,data.frame,attach,detach,describe,datadensity,plot.data.frame,hist.data.frame,naclus,factor,label,units,names,expand.grid,summary.formula,summary.data.frame,casefold,edit,page,plot.data.frame,Cs,combine.levels,upData

Examples

## Not run: # First, we do steps that create or manipulate the data# frame in its entirety.  For S-Plus, these are done with# .Data in search position one (the default at the# start of the session).## -----------------------------------------------------------------------# Step 1: Create initial draft of data frame# # We usually begin by importing a dataset from# # another application.  ASCII files may be imported# using the scan and read.table functions.  SAS# datasets may be imported using the Hmisc sas.get# function (which will carry more attributes from# SAS than using File \dots  Import) from the GUI# menus.  But for most applications (especially# Excel), File \dots Import will suffice.  If using# the GUI, it is often best to provide variable# names during the import process, using the Options# tab, rather than renaming all fields later Of# course, if the data to be imported already have# field names (e.g., in Excel), let S use those# automatically.  If using S-Plus, you can use a# command to execute File \dots  Import, e.g.:import.data(FileName = "/windows/temp/fev.asc",            FileType = "ASCII", DataFrame = "FEV")# Here we name the new data frame FEV rather than# fev, because we wanted to distinguish a variable# in the data frame named fev from the data frame# name.  For S-Plus the command will look# instead like the following:FEV <- importData("/tmp/fev.asc")# -----------------------------------------------------------------------# Step 2: Clean up data frame / make it be more# efficiently stored# # Unless using sas.get to import your dataset# (sas.get already stores data efficiently), it is# usually a good idea to run the data frame through# the Hmisc cleanup.import function to change# numeric variables that are always whole numbers to# be stored as integers, the remaining numerics to# single precision, strange values from Excel to# NAs, and character variables that always contain# legal numeric values to numeric variables.# cleanup.import typically halves the size of the# data frame.  If you do not specify any parameters# to cleanup.import, the function assumes that no# numeric variable needs more than 7 significant# digits of precision, so all non-integer-valued# variables will be converted to single precision.FEV <- cleanup.import(FEV)# -----------------------------------------------------------------------# Step 3: Make global changes to the data frame# # A data frame has attributes that are "external" to# its variables.  There are the vector of its# variable names ("names" attribute), the# observation identifiers ("row.names"), and the# "class" (whose value is "data.frame").  The# "names" attribute is the one most commonly in need# of modification.  If we had wanted to change all# the variable names to lower case, we could have# specified lowernames=TRUE to the cleanup.import# invocation above, or typenames(FEV) <- casefold(names(FEV))# The upData function can also be used to change# variable names in two ways (see below).# To change names in a non-systematic way we use# other options.  Under Windows/NT the most# straigtforward approach is to change the names# interactively.  Click on the data frame in the# left panel of the Object Browser, then in the# right pane click twice (slowly) on a variable.# Use the left arrow and other keys to edit the# name.  Click outside that name field to commit the# change.  You can also rename columns while in a# Data Sheet.  To instead use programming commands# to change names, use something like:names(FEV)[6] <- 'smoke'   # assumes you know the positions!  names(FEV)[names(FEV)=='smoking'] <- 'smoke' names(FEV) <- edit(names(FEV))# The last example is useful if you are changing# many names.  But none of the interactive# approaches such as edit() are handy if you will be# re-importing the dataset after it is updated in# its original application.  This problem can be# addressed by saving the new names in a permanent# vector in .Data:new.names <- names(FEV)# Then if the data are re-imported, you can typenames(FEV) <- new.names# to rename the variables.# -----------------------------------------------------------------------# Step 4: Delete unneeded variables# # To delete some of the variables, you can# right-click on variable names in the Object# Browser's right pane, then select Delete.  You can# also set variables to have NULL values, which# causes the system to delete them.  We don't need# to delete any variables from FEV but suppose we# did need to delete some from mydframe.mydframe$x1 <- NULL mydframe$x2 <- NULLmydframe[c('age','sex')] <- NULL   # delete 2 variables mydframe[Cs(age,sex)]    <- NULL   # same thing# The last example uses the Hmisc short-cut quoting# function Cs.  See also the drop parameter to upData.# -----------------------------------------------------------------------# Step 5: Make changes to individual variables#         within the data frame# # After importing data, the resulting variables are# seldom self - documenting, so we commonly need to# change or enhance attributes of individual# variables within the data frame.# # If you are only changing a few variables, it is# efficient to change them directly without# attaching the entire data frame.FEV$sex   <- factor(FEV$sex,   0:1, c('female','male')) FEV$smoke <- factor(FEV$smoke, 0:1,                     c('non-current smoker','current smoker')) units(FEV$age)    <- 'years'units(FEV$fev)    <- 'L' label(FEV$fev)    <- 'Forced Expiratory Volume' units(FEV$height) <- 'inches'# When changing more than one or two variables it is# more convenient change the data frame using the# Hmisc upData function.FEV2 <- upData(FEV,  rename=c(smoking='smoke'),   # omit if renamed above  drop=c('var1','var2'),  levels=list(sex  =list(female=0,male=1),              smoke=list('non-current smoker'=0,                         'current smoker'=1)),  units=list(age='years', fev='L', height='inches'),  labels=list(fev='Forced Expiratory Volume'))# An alternative to levels=list(\dots) is for example# upData(FEV, sex=factor(sex,0:1,c('female','male'))).# # Note that we saved the changed data frame into a# new data frame FEV2.  If we were confident of the# correctness of our changes we could have stored# the new data frame on top of the old one, under# the original name FEV.# -----------------------------------------------------------------------# Step 6:  Check the data frame# # The Hmisc describe function is perhaps the first# function that should be used on the new data# frame.  It provides documentation of all the# variables and the frequency tabulation, counts of# NAs,  and 5 largest and smallest values are# helpful in detecting data errors.  Typing# describe(FEV) will write the results to the# current output window.  To put the results in a# new window that can persist, even upon exiting# S, we use the page function.  The describe# output can be minimized to an icon but kept ready# for guiding later steps of the analysis.page(describe(FEV2), multi=TRUE) # multi=TRUE allows that window to persist while# control is returned to other windows# The new data frame is OK.  Store it on top of the# old FEV and then use the graphical user interface# to delete FEV2 (click on it and hit the Delete# key) or type rm(FEV2) after the next statement.FEV <- FEV2# Next, we can use a variety of other functions to# check and describe all of the variables.  As we# are analyzing all or almost all of the variables,# this is best done without attaching the data# frame.  Note that plot.data.frame plots inverted# CDFs for continuous variables and dot plots# showing frequency distributions of categorical# ones.summary(FEV)# basic summary function (summary.data.frame) plot(FEV)                # plot.data.frame datadensity(FEV)         # rug plots and freq. bar charts for all var.hist.data.frame(FEV)     # for variables having > 2 values by(FEV, FEV$smoke, summary)  # use basic summary function with stratification# -----------------------------------------------------------------------# Step 7:  Do detailed analyses involving individual#          variables# # Analyses based on the formula language can use# data= so attaching the data frame may not be# required.  This saves memory.  Here we use the# Hmisc summary.formula function to compute 5# statistics on height, stratified separately by age# quartile and by sex.options(width=80) summary(height ~ age + sex, data=FEV,        fun=function(y)c(smean.sd(y),                         smedian.hilow(y,conf.int=.5)))# This computes mean height, S.D., median, outer quartilesfit <- lm(height ~ age*sex, data=FEV) summary(fit)# For this analysis we could also have attached the# data frame in search position 2.  For other# analyses, it is mandatory to attach the data frame# unless FEV$ prefixes each variable name.# Important: DO NOT USE attach(FEV, 1) or# attach(FEV, pos=1, \dots) if you are only analyzing# and not changing the variables, unless you really# need to avoid conflicts with variables in search# position 1 that have the same names as the# variables in FEV.  Attaching into search position# 1 will cause S-Plus to be more of a memory hog.attach(FEV)# Use e.g. attach(FEV[,Cs(age,sex)]) if you only# want to analyze a small subset of the variables# Use e.g. attach(FEV[FEV$sex=='male',]) to# analyze a subset of the observationssummary(height ~ age + sex,        fun=function(y)c(smean.sd(y),          smedian.hilow(y,conf.int=.5)))fit <- lm(height ~ age*sex)# Run generic summary function on height and fev, # stratified by sexby(data.frame(height,fev), sex, summary)# Cross-classify into 4 sex x smoke groupsby(FEV, list(sex,smoke), summary)# Plot 5 quantiless <- summary(fev ~ age + sex + height,              fun=function(y)quantile(y,c(.1,.25,.5,.75,.9)))plot(s, which=1:5, pch=c(1,2,15,2,1), #pch=c('=','[','o',']','='),      main='A Discovery', xlab='FEV')# Use the nonparametric bootstrap to compute a # 0.95 confidence interval for the population mean fevsmean.cl.boot(fev)    # in Hmisc# Use the Statistics \dots Compare Samples \dots One Sample # keys to get a normal-theory-based C.I.  Then do it # more manually.  The following method assumes that # there are no NAs in fevsd <- sqrt(var(fev))xbar <- mean(fev)xbarsdn <- length(fev)qt(.975,n-1)     # prints 0.975 critical value of t dist. with n-1 d.f.xbar + c(-1,1)*sd/sqrt(n)*qt(.975,n-1)   # prints confidence limits# Fit a linear model# fit <- lm(fev ~ other variables \dots)detach()# The last command is only needed if you want to# start operating on another data frame and you want# to get FEV out of the way.# -----------------------------------------------------------------------# Creating data frames from scratch# # Data frames can be created from within S.  To# create a small data frame containing ordinary# data, you can use something likedframe <- data.frame(age=c(10,20,30),                      sex=c('male','female','male'),                     stringsAsFactors=TRUE)# You can also create a data frame using the Data# Sheet.  Create an empty data frame with the# correct variable names and types, then edit in the# data.dd <- data.frame(age=numeric(0),sex=character(0),                 stringsAsFactors=TRUE)# The sex variable will be stored as a factor, and# levels will be automatically added to it as you# define new values for sex in the Data Sheet's sex# column.# # When the data frame you need to create is defined# by systematically varying variables (e.g., all# possible combinations of values of each variable),# the expand.grid function is useful for quickly# creating the data.  Then you can add# non-systematically-varying variables to the object# created by expand.grid, using programming# statements or editing the Data Sheet.  This# process is useful for creating a data frame# representing all the values in a printed table.# In what follows we create a data frame# representing the combinations of values from an 8# x 2 x 2 x 2 (event x method x sex x what) table,# and add a non-systematic variable percent to the# data.jcetable <- expand.grid( event=c('Wheezing at any time',         'Wheezing and breathless',         'Wheezing without a cold',         'Waking with tightness in the chest',         'Waking with shortness of breath',         'Waking with an attack of cough',         'Attack of asthma',         'Use of medication'), method=c('Mail','Telephone'),  sex=c('Male','Female'), what=c('Sensitivity','Specificity'))jcetable$percent <- c(756,618,706,422,356,578,289,333,  576,421,789,273,273,212,212,212,  613,763,713,403,377,541,290,226,  613,684,632,290,387,613,258,129,  656,597,438,780,732,679,938,919,  714,600,494,877,850,703,963,987,  755,420,480,794,779,647,956,941,  766,423,500,833,833,604,955,986) / 10# In jcetable, event varies most rapidly, then# method, then sex, and what.## End(Not run)

Representativeness of Observations in a Data Set

Description

These functions are intended to be used to describe how well a givenset of new observations (e.g., new subjects) were represented in adataset used to develop a predictive model.ThedataRep function forms a data frame that contains all the uniquecombinations of variable values that existed in a given set ofvariable values. Cross–classifications of values are created usingexact values of variables, so for continuous numeric variables it isoften necessary to round them to the nearestv and to possiblycurtail the values to some lower and upper limit before rounding.Herev denotes a numeric constant specifying the matching tolerancethat will be used.dataRep also stores marginal distributionsummaries for all the variables. For numeric variables, all 101percentiles are stored, and for all variables, the frequencydistributions are also stored (frequencies are computed after anyrounding and curtailment of numeric variables). For the purposes ofrounding and curtailing, theroundN function is provided. Aprintmethod will summarize the calculations made bydataRep, and iflong=TRUE all unique combinations of values and their frequencies inthe original dataset are printed.

Thepredict method fordataRep takes a new data frame havingvariables named the same as the original ones (but whose factor levelsare not necessarily in the same order) and examines the collapsedcross-classifications created bydataRep to find how manyobservations were similar to each of the new observations after anyrounding or curtailment of limits is done.predict also does somecalculations to describe how the variable values of the newobservations "stack up" against the marginal distributions of theoriginal data. For categorical variables, the percent of observationshaving a given variable with the value of the new observation (afterrounding for variables that were throughroundN in the formula giventodataRep) is computed. For numeric variables, the percentile ofthe original distribution in which the current value falls will becomputed. For this purpose, the data are not rounded because the 101original percentiles were retained; linear interpolation is used toestimate percentiles for values between two tabulated percentiles.The lowest marginal frequency of matching values across all variablesis also computed. For example, if an age, sex combination matches 10subjects in the original dataset but the age value matches 100 ages(after rounding) and the sex value matches the sex code of 300observations, the lowest marginal frequency is 100, which is a "bestcase" upper limit for multivariable matching. I.e., matching on allvariables has to result on a lower frequency than this amount.Aprint method for the output ofpredict.dataRep prints allcalculations done bypredict by default. Calculations can beselectively suppressed.

Usage

dataRep(formula, data, subset, na.action)roundN(x, tol=1, clip=NULL)## S3 method for class 'dataRep'print(x, long=FALSE, ...)## S3 method for class 'dataRep'predict(object, newdata, ...)## S3 method for class 'predict.dataRep'print(x, prdata=TRUE, prpct=TRUE, ...)

Arguments

Value

dataRep returns a list of class"dataRep" containing the collapseddata frame and frequency counts along with marginal distributioninformation.predict returns an object of class"predict.dataRep"containing information determined by matching observations innewdata with the original (collapsed) data.

Side Effects

print.dataRep prints.

Author(s)

Frank Harrell
Department of Biostatistics
Vanderbilt University School of Medicine
fh@fharrell.com

See Also

round,table

Examples

set.seed(13)num.symptoms <- sample(1:4, 1000,TRUE)sex <- factor(sample(c('female','male'), 1000,TRUE))x    <- runif(1000)x[1] <- NAtable(num.symptoms, sex, .25*round(x/.25))d <- dataRep(~ num.symptoms + sex + roundN(x,.25))print(d, long=TRUE)predict(d, data.frame(num.symptoms=1:3, sex=c('male','male','female'),                      x=c(.03,.5,1.5)))

Design Effect and Intra-cluster Correlation

Description

Computes the Kish design effect and corresponding intra-cluster correlationfor a single cluster-sampled variable

Usage

deff(y, cluster)

Arguments

Value

a vector with named elementsn (total number of non-missingobservations),clusters (number of clusters after deletingmissing data),rho(intra-cluster correlation), anddeff(design effect).

Author(s)

Frank Harrell
Department of Biostatistics
Vanderbilt University
fh@fharrell.com

See Also

bootcov,robcov

Examples

set.seed(1)blood.pressure <- rnorm(1000, 120, 15)clinic <- sample(letters, 1000, replace=TRUE)deff(blood.pressure, clinic)

Concise Statistical Description of a Vector, Matrix, Data Frame,or Formula

Description

describe is a generic method that invokesdescribe.data.frame,describe.matrix,describe.vector, ordescribe.formula.describe.vector is the basic function for handling a single variable.This function determines whether the variable is character, factor,category, binary, discrete numeric, and continuous numeric, and printsa concise statistical summary according to each. A numeric variable isdeemed discrete if it has <= 10 distinct values. In this case,quantiles are not printed. A frequency table is printed for any non-binary variable if it has no more than 20 distinctvalues. For any variable for which the frequency table is not printed,the 5 lowest and highest values are printed. This behavior can beoverriden for long character variables with many levels using thelistunique parameter, to get a complete tabulation.

describe is especially useful fordescribing data frames created by*.get, as labels, formats,value labels, and (in the case ofsas.get) frequencies of specialmissing values are printed.

For a binary variable, the sum (number of 1's) and mean (proportion of1's) are printed. If the first argument is a formula, a model frameis created and passed to describe.data.frame. If a variableis of class"impute", a count of the number of imputed values isprinted. If a date variable has an attributepartial.date(this is set up bysas.get), counts of how many partial dates areactually present (missing month, missing day, missing both) are also presented.If a variable was created by the special-purpose functionsubsti (whichsubstitutes values of a second variable if the first variable is NA),the frequency table of substitutions is also printed.

For numeric variables,describe adds an item calledInfowhich is a relative information measure using the relative efficiency ofa proportional odds/Wilcoxon test on the variable relative to the sametest on a variable that has no ties.Info is related to howcontinuous the variable is, and ties are less harmful the more untiedvalues there are. The formula forInfo is one minus the sum ofthe cubes of relative frequencies of values divided by one minus thesquare of the reciprocal of the sample size. The lowest informationcomes from a variable having only one distinct value following by ahighly skewed binary variable.Info is reported totwo decimal places.

A latex method exists for converting thedescribe object to aLaTeX file. For numeric variables having more than 20 distinct values,describe saves in its returned object the frequencies of 100evenly spaced bins running from minimum observed value to the maximum.When there are less than or equal to 20 distinct values, the originalvalues are maintained.latex andhtml insert a spike histogram displaying thesefrequency counts in the tabular material using the LaTeX pictureenvironment. For example output seehttps://hbiostat.org/doc/rms/book/chapter7edition1.pdf.Note that the latex method assumes you have the following stylesinstalled in your latex installation: setspace and relsize.

Thehtml method mimics the LaTeX output. This is useful in thecontext of Quarto/Rmarkdown html and html notebook output.Ifoptions(prType='html') is in effect, callingprint onan object that is the result of runningdescribe on a data framewill result in rendering the HTML version. If run from the console abrowser window will open. Whenwhich is specified toprint, whether or notprType='html' is in effect, agt package html table will be produced containing only the types of variables requested. Whenwhich='both' a list withelement namesContinuous andCategorical is produced,making it convenient for the user to print as desired, or to pass thelist directed to theqreportmaketabs function when using Quarto.

Theplot method is fordescribe objects run on dataframes. It produces spike histograms for a graphic ofcontinuous variables and a dot chart for categorical variables, showingcategory proportions. The graphic format isggplot2 if the userhas not setoptions(grType='plotly') or has set thegrTypeoption to something other than'plotly'. Otherwiseplotlygraphics that are interactive are produced, and these can be placed intoan Rmarkdown html notebook. The user must install theplotlypackage for this to work. When the use hovers the mouse over a bin fora raw data value, the actual value will pop-up (formatted usingdigits). When the user hovers over the minimum data value, mostof the information calculated bydescribe will pop up. For eachvariable, the number of missing values is used to assign the color tothe histogram or dot chart, and a legend is drawn. Color is not used ifthere are no missing values in any variable. For categorical variables,hovering over the leftmost point for a variable displays details, andfor all points proportions, numerators, and denominators are displayedin the popup. If both continuous and categorical variables are presentandwhich='both' is specified, theplot method returns anunclassedlist containing two objects, named'Categorical'and'Continuous', in that order.

Sample weights may be specified to any of the functions, resultingin weighted means, quantiles, and frequency tables.

Note: As discussed in Cox and Longton (2008), Stata Technical Bulletin 8(4)pp. 557, the term "unique" has been replaced with "distinct" in theoutput (but not in parameter names).

Whenweights are not used, the pseudomedian and Gini's mean difference are computed fornumeric variables. The pseudomedian is labeledpMedian and is the median of all possible pairwise averages. It is a robust and efficient measure of location that equals the mean and median for symmetric distributions. It is also called the Hodges-Lehmann one-sample estimator. Gini's mean difference is a robust measure of dispersion that is themean absolute difference between any pairs of observations. In simpleoutput Gini's difference is labeledGmd.

formatdescribeSingle is a service function forlatex,html, andprint methods for single variables that is notintended to be called by the user.

Usage

## S3 method for class 'vector'describe(x, descript, exclude.missing=TRUE, digits=4,         listunique=0, listnchar=12,         weights=NULL, normwt=FALSE, minlength=NULL, shortmChoice=TRUE,         rmhtml=FALSE, trans=NULL, lumptails=0.01, ...)## S3 method for class 'matrix'describe(x, descript, exclude.missing=TRUE, digits=4, ...)## S3 method for class 'data.frame'describe(x, descript, exclude.missing=TRUE,    digits=4, trans=NULL, ...)## S3 method for class 'formula'describe(x, descript, data, subset, na.action,    digits=4, weights, ...)## S3 method for class 'describe'print(x, which = c('both', 'categorical', 'continuous'), ...)## S3 method for class 'describe'latex(object, title=NULL,      file=paste('describe',first.word(expr=attr(object,'descript')),'tex',sep='.'),      append=FALSE, size='small', tabular=TRUE, greek=TRUE,      spacing=0.7, lspace=c(0,0), ...)## S3 method for class 'describe.single'latex(object, title=NULL, vname,      file, append=FALSE, size='small', tabular=TRUE, greek=TRUE,      lspace=c(0,0), ...)## S3 method for class 'describe'html(object, size=85, tabular=TRUE,      greek=TRUE, scroll=FALSE, rows=25, cols=100, ...)## S3 method for class 'describe.single'html(object, size=85,      tabular=TRUE, greek=TRUE, ...)formatdescribeSingle(x, condense=c('extremes', 'frequencies', 'both', 'none'),           lang=c('plain', 'latex', 'html'), verb=0, lspace=c(0, 0),           size=85, ...)## S3 method for class 'describe'plot(x, which=c('both', 'continuous', 'categorical'),                          what=NULL,                          sort=c('ascending', 'descending', 'none'),                          n.unique=10, digits=5, bvspace=2, ...)

Arguments

Details

Ifoptions(na.detail.response=TRUE)has been set andna.action is"na.delete" or"na.keep", summary statistics onthe response variable are printed separately for missing and non-missingvalues of each predictor. The default summary function returnsthe number of non-missing response values and the mean of the lastcolumn of the response values, with anames attribute ofc("N","Mean"). When the response is aSurv object and the mean is used, this willresult in the crude proportion of events being used to summarizethe response. The actual summary function can be designated throughoptions(na.fun.response = "function name").

If you are modifying LaTexparskip or certain other parameters,you may need to shrink the area aroundtabular andverbatim environments produced bylatex.describe. You cando this using for example\usepackage{etoolbox}\makeatletter\preto{\@verbatim}{\topsep=-1.4pt\partopsep=0pt}\preto{\@tabular}{\parskip=2pt\parsep=0pt}\makeatother in the LaTeX preamble.

Multiple choice (mChoice) variables'describe output renders well in html but not when included in aQuarto document.

Value

a list containing elementsdescript,counts,values. The list is of classdescribe. If the inputobject was a matrix or a data frame, the list is a list of lists, one list for each variableanalyzed.latex returns a standardlatex object. For numericvariables having at least 20 distinct values, an additional componentintervalFreq. This component is a list with two elements,range(containing two values) andcount, a vector of 100 integer frequencycounts.print withwhich= returns a 'gt' table object.The user can modify the table by piping formatting changes, columnremovals, and other operations, before final rendering.

Author(s)

Frank Harrell
Vanderbilt University
fh@fharrell.com

See Also

spikecomp,sas.get,quantile,GiniMd,pMedian,table,summary,model.frame.default,naprint,lapply,tapply,Surv,na.delete,na.keep,na.detail.response,latex

Examples

set.seed(1)describe(runif(200),dig=2)    #single variable, continuous                              #get quantiles .05,.10,\dotsdfr <- data.frame(x=rnorm(400),y=sample(c('male','female'),400,TRUE))describe(dfr)## Not run: options(grType='plotly')d <- describe(mydata)p <- plot(d)   # create plots for both types of variablesp[[1]]; p[[2]] # or p$Categorical; p$Continuousplotly::subplot(p[[1]], p[[2]], nrows=2)  # plot both in oneplot(d, which='categorical')    # categorical onesd <- sas.get(".","mydata",special.miss=TRUE,recode=TRUE)describe(d)      #describe entire data frameattach(d, 1)describe(relig)  #Has special missing values .D .F .M .R .T                 #attr(relig,"label") is "Religious preference"#relig : Religious preference  Format:relig#    n missing  D  F M R T distinct # 4038     263 45 33 7 2 1        8##0:none (251, 6%), 1:Jewish (372, 9%), 2:Catholic (1230, 30%) #3:Jehovah's Witnes (25, 1%), 4:Christ Scientist (7, 0%) #5:Seventh Day Adv (17, 0%), 6:Protestant (2025, 50%), 7:other (111, 3%) # Method for describing part of a data frame: describe(death.time ~ age*sex + rcs(blood.pressure)) describe(~ age+sex) describe(~ age+sex, weights=freqs)  # weighted analysis fit <- lrm(y ~ age*sex + log(height)) describe(formula(fit)) describe(y ~ age*sex, na.action=na.delete)   # report on number deleted for each variable options(na.detail.response=TRUE)  # keep missings separately for each x, report on dist of y by x=NA describe(y ~ age*sex) options(na.fun.response="quantile") describe(y ~ age*sex)   # same but use quantiles of y by x=NA d <- describe(my.data.frame) d$age                   # print description for just age d[c('age','sex')]       # print description for two variables d[sort(names(d))]       # print in alphabetic order by var. names d2 <- d[20:30]          # keep variables 20-30 page(d2)                # pop-up window for these variables# Test date/time formats and suppression of times when they don't vary library(chron) d <- data.frame(a=chron((1:20)+.1),                 b=chron((1:20)+(1:20)/100),                 d=ISOdatetime(year=rep(2003,20),month=rep(4,20),day=1:20,                               hour=rep(11,20),min=rep(17,20),sec=rep(11,20)),                 f=ISOdatetime(year=rep(2003,20),month=rep(4,20),day=1:20,                               hour=1:20,min=1:20,sec=1:20),                 g=ISOdate(year=2001:2020,month=rep(3,20),day=1:20)) describe(d)# Make a function to run describe, latex.describe, and use the kdvi# previewer in Linux to view the result and easily make a pdf file ldesc <- function(data) {  options(xdvicmd='kdvi')  d <- describe(data, desc=deparse(substitute(data)))  dvi(latex(d, file='/tmp/z.tex'), nomargins=FALSE, width=8.5, height=11) } ldesc(d)## End(Not run)

Discrete Vector tools

Description

discrete creates a discrete vector which is distinct from acontinuous vector, or a factor/ordered vector.The other function are tools for manipulating descrete vectors.

Usage

as.discrete(x, ...)## Default S3 method:as.discrete(x, ...)discrete(x, levels = sort(unique.default(x), na.last = TRUE), exclude = NA)## S3 replacement method for class 'discrete'x[...] <- value## S3 method for class 'discrete'x[..., drop = FALSE]## S3 method for class 'discrete'x[[i]]is.discrete(x)## S3 replacement method for class 'discrete'is.na(x) <- value## S3 replacement method for class 'discrete'length(x) <- value

Arguments

Details

as.discrete converts a vector into a discrete vector.

discrete creates a discrete vector from provided values.

is.discrete tests to see if the vector is a discrete vector.

Value

as.discrete,discrete returns a vector ofdiscrete type.

is.discrete returan logicalTRUE if the vector is ofclass discrete other wise it returnsFALSE.

Author(s)

Charles Dupont

See Also

[[,[,factor

Examples

a <- discrete(1:25)ais.discrete(a)b <- as.discrete(2:4)b

Enhanced Dot Chart

Description

dotchart2 is an enhanced version of thedotchart function with several new options.

Usage

dotchart2(data, labels, groups=NULL, gdata=NA, horizontal=TRUE, pch=16,          xlab='', ylab='', xlim=NULL, auxdata, auxgdata=NULL, auxtitle,          lty=1, lines=TRUE, dotsize = .8,          cex = par("cex"), cex.labels = cex,          cex.group.labels = cex.labels*1.25, sort.=TRUE,       add=FALSE, dotfont=par('font'), groupfont=2,       reset.par=add, xaxis=TRUE, width.factor=1.1,          lcolor='gray', leavepar=FALSE,          axisat=NULL, axislabels=NULL, ...)

Arguments

Side Effects

dotchart will leavepar altered ifreset.par=FALSE.

Author(s)

Frank Harrell
Department of Biostatistics
Vanderbilt University
fh@fharrell.com

See Also

dotchart

Examples

set.seed(135)maj <- factor(c(rep('North',13),rep('South',13)))g <- paste('Category',rep(letters[1:13],2))n <- sample(1:15000, 26, replace=TRUE)y1 <- runif(26)y2 <- pmax(0, y1 - runif(26, 0, .1))dotchart2(y1, g, groups=maj, auxdata=n, auxtitle='n', xlab='Y')dotchart2(y2, g, groups=maj, pch=17, add=TRUE)## Compare with dotchart function (no superpositioning or auxdata allowed):## dotchart(y1, g, groups=maj, xlab='Y')## To plot using a transformed scale add for example## axisat=sqrt(pretty(y)), axislabels=pretty(y)

Enhanced Version of dotchart Function

Description

These are adaptations of the R dotchart function that sorts categoriestop to bottom, addsauxdata andauxtitle arguments to putextra information in the right margin, and fordotchart3 addsargumentscex.labels,cex.group.labels, andgroupfont. By default, group headings are in a larger, boldfont.dotchart3 also cuts a bit of white space from the top andbottom of the chart. The most significant change, however, is in howx is interpreted. Columns ofx no longer provide analternate way to define groups. Instead, they define superpositionedvalues. This is useful for showing three quartiles, for example. Goingalong with this change, fordotchart3pch can now be avector specifying symbols to use going across columns ofx.x was changed in this way because to put multiple points on aline (e.g., quartiles) and keeping track ofpar() parameters whendotchart2 was called withadd=TRUE was cumbersome.dotchart3 changes the margins to account for horizontal labels.

dotchartp is a version ofdotchart3 for making the chartwith theplotly package.

summaryD creates aggregate data usingsummarize andcallsdotchart3 with suitable arguments to summarize data bymajor and minor categories. Ifoptions(grType='plotly') is ineffect and theplotly package is installed,summaryD usesdotchartp instead ofdotchart3.

summaryDp is a streamlinedsummaryD-like function thatuses thedotchartpl function to render aplotly graphic.It is used to compute summary statistics stratified separately by aseries of variables.

Usage

dotchart3(x, labels = NULL, groups = NULL, gdata = NULL,          cex = par("cex"), pch = 21, gpch = pch, bg = par("bg"),          color = par("fg"), gcolor = par("fg"), lcolor = "gray",          xlim = range(c(x, gdata), na.rm=TRUE), main = NULL, xlab = NULL,          ylab = NULL, auxdata = NULL, auxtitle = NULL, auxgdata=NULL,          axisat=NULL, axislabels=NULL,          cex.labels = cex, cex.group.labels = cex.labels * 1.25,          cex.auxdata=cex, groupfont = 2,          auxwhere=NULL, height=NULL, width=NULL, ...)dotchartp(x, labels = NULL, groups = NULL, gdata = NULL,            xlim = range(c(x, gdata), na.rm=TRUE), main=NULL,            xlab = NULL, ylab = '', auxdata=NULL, auxtitle=NULL,            auxgdata=NULL, auxwhere=c('right', 'hover'),            symbol='circle', col=colorspace::rainbow_hcl,            legendgroup=NULL,            axisat=NULL, axislabels=NULL, sort=TRUE, digits=4, dec=NULL,            height=NULL, width=700, layoutattr=FALSE, showlegend=TRUE, ...) summaryD(formula, data=NULL, fun=mean, funm=fun,         groupsummary=TRUE, auxvar=NULL, auxtitle='',         auxwhere=c('hover', 'right'),         vals=length(auxvar) > 0, fmtvals=format,         symbol=if(use.plotly) 'circle' else 21,         col=if(use.plotly) colorspace::rainbow_hcl else 1:10,         legendgroup=NULL,         cex.auxdata=.7, xlab=v[1], ylab=NULL,         gridevery=NULL, gridcol=gray(.95), sort=TRUE, ...)summaryDp(formula,          fun=function(x) c(Mean=mean(x, na.rm=TRUE),                            N=sum(! is.na(x))),          overall=TRUE, xlim=NULL, xlab=NULL,          data=NULL, subset=NULL, na.action=na.retain,          ncharsmax=c(50, 30),          digits=4, ...)

Arguments

Value

the function returns invisibly

Author(s)

Frank Harrell

See Also

dotchart,dotchart2,summarize,rlegend

Examples

set.seed(135)maj <- factor(c(rep('North',13),rep('South',13)))g <- paste('Category',rep(letters[1:13],2))n <- sample(1:15000, 26, replace=TRUE)y1 <- runif(26)y2 <- pmax(0, y1 - runif(26, 0, .1))dotchart3(cbind(y1,y2), g, groups=maj, auxdata=n, auxtitle='n',          xlab='Y', pch=c(1,17))## Compare with dotchart function (no superpositioning or auxdata allowed):## dotchart(y1, g, groups=maj, xlab='Y')## Not run: dotchartp(cbind(y1, y2), g, groups=maj, auxdata=n, auxtitle='n',          xlab='Y', gdata=cbind(c(0,.1), c(.23,.44)), auxgdata=c(-1,-2),          symbol=c('circle', 'line-ns-open'))summaryDp(sbp ~ region + sex + race + cut2(age, g=5), data=mydata)## End(Not run)## Put options(grType='plotly') to have the following use dotchartp## (rlegend will not apply)## Add argument auxwhere='hover' to summaryD or dotchartp to put## aux info in hover text instead of right marginsummaryD(y1 ~ maj + g, xlab='Mean')summaryD(y1 ~ maj + g, groupsummary=FALSE)summaryD(y1 ~ g, fmtvals=function(x) sprintf('%4.2f', x))Y <- cbind(y1, y2)   # summaryD cannot handle cbind(...) ~ ...summaryD(Y  ~ maj + g, fun=function(y) y[1,], symbol=c(1,17))rlegend(.1, 26, c('y1','y2'), pch=c(1,17))summaryD(y1 ~ maj, fun=function(y) c(Mean=mean(y), n=length(y)),         auxvar='n', auxtitle='N')

Enhanced Version of dotchart Function for plotly

Description

This function produces aplotly interactive graphic and acceptsa different format of data input than the otherdotchartfunctions. It was written to handle a hierarchical data structureincluding strata that further subdivide the main classes. Strata,indicated by themult variable, are shown on the samehorizontal line, and if the variablebig isFALSE willappear slightly below the main line, using smaller symbols, and havingsome transparency. This is intended to handle output such as thatfrom thesummaryP function when there is a superpositioningvariablegroup and a stratification variablemult,especially when the data have been run through theaddMarginalfunction to createmult categories labelled"All" forwhich the user will specifybig=TRUE to indicate non-stratifiedestimates (stratified only ongroup) to emphasize.

When viewing graphics that usedmult andbig, the usercan click on the legends for the small points forgroups tovanish the finely stratified estimates.

Whengroup is used bymult andbig are not, andwhen thegroup variable has exactly two distinct values, youcan specifyrefgroup to get the difference between twoproportions in addition to the individual proportions. The individualproportions are plotted, but confidence intervals for the differenceare shown in hover text and half-width confidence intervals for thedifference, centered at the midpoint of the proportions, are shown.These have the property of intersecting the two proportions if andonly if there is no significant difference at the1 - conf.intlevel.

Specifyfun=exp andifun=log if estimates and confidencelimits are on the log scale. Make sure that zeros were prevented inthe original calculations. For exponential hazard rates this can beaccomplished by replacing event counts of 0 with 0.5.

Usage

dotchartpl(x, major=NULL, minor=NULL, group=NULL, mult=NULL,           big=NULL, htext=NULL, num=NULL, denom=NULL,           numlabel='', denomlabel='',           fun=function(x) x, ifun=function(x) x, op='-',           lower=NULL, upper=NULL,           refgroup=NULL, sortdiff=TRUE, conf.int=0.95,           minkeep=NULL, xlim=NULL, xlab='Proportion',           tracename=NULL, limitstracename='Limits',           nonbigtracename='Stratified Estimates',           dec=3, width=800, height=NULL,           col=colorspace::rainbow_hcl)

Arguments

Value

aplotly object. An attributelevelsRemoved isadded ifminkeep is used and any categories were omitted fromthe plot as a result. This is a character vector with categoriesremoved. Ifmajor is present, the strings are of the formmajor:minor

Author(s)

Frank Harrell

See Also

dotchartp

Examples

## Not run: set.seed(1)d <- expand.grid(major=c('Alabama', 'Alaska', 'Arkansas'),                 minor=c('East', 'West'),                 group=c('Female', 'Male'),                 city=0:2)n <- nrow(d)d$num   <- round(100*runif(n))d$denom <- d$num + round(100*runif(n))d$x     <- d$num / d$denomd$lower <- d$x - runif(n)d$upper <- d$x + runif(n)with(d, dotchartpl(x, major, minor, group, city, lower=lower, upper=upper,            big=city==0, num=num, denom=denom, xlab='x'))# Show half-width confidence intervals for Female - Male differences# after subsetting the data to have only one record per# state/region/groupd <- subset(d, city == 0)with(d, dotchartpl(x, major, minor, group, num=num, denom=denom,            lower=lower, upper=upper, refgroup='Male'))n <- 500set.seed(1)d <- data.frame(  race         = sample(c('Asian', 'Black/AA', 'White'), n, TRUE),  sex          = sample(c('Female', 'Male'), n, TRUE),  treat        = sample(c('A', 'B'), n, TRUE),  smoking      = sample(c('Smoker', 'Non-smoker'), n, TRUE),  hypertension = sample(c('Hypertensive', 'Non-Hypertensive'), n, TRUE),  region       = sample(c('North America','Europe','South America',                          'Europe', 'Asia', 'Central America'), n, TRUE))d <- upData(d, labels=c(race='Race', sex='Sex'))dm <- addMarginal(d, region)s <- summaryP(race + sex + smoking + hypertension ~                region + treat,  data=dm)s$region <- ifelse(s$region == 'All', 'All Regions', as.character(s$region))with(s,  dotchartpl(freq / denom, major=var, minor=val, group=treat, mult=region,            big=region == 'All Regions', num=freq, denom=denom))s2 <- s[- attr(s, 'rows.to.exclude1'), ]with(s2,      dotchartpl(freq / denom, major=var, minor=val, group=treat, mult=region,                big=region == 'All Regions', num=freq, denom=denom))# Note these plots can be created by plot.summaryP when options(grType='plotly')# Plot hazard rates and ratios with confidence limits, on log scaled <- data.frame(tx=c('a', 'a', 'b', 'b'),                event=c('MI', 'stroke', 'MI', 'stroke'),                count=c(10, 5, 5, 2),                exposure=c(1000, 1000, 900, 900))# There were no zero event counts in this dataset.  In general we# want to handle that, hence the 0.5 belowd <- upData(d, hazard = pmax(0.5, count) / exposure,               selog  = sqrt(1. / pmax(0.5, count)),               lower  = log(hazard) - 1.96 * selog,               upper  = log(hazard) + 1.96 * selog)with(d,     dotchartpl(log(hazard), minor=event, group=tx, num=count, denom=exposure,                lower=lower, upper=upper,                fun=exp, ifun=log, op='/',                numlabel='events', denomlabel='years',                refgroup='a', xlab='Events Per Person-Year'))## End(Not run)

Dual Standard Deviations

Description

Computes one standard deviation for the lower half of the distribution of a numeric vector and another SD for the upper half. By default the center of the distribution for purposes of splitting into "halves" is the mean. The user may override this withcenter. When splitting into halves, observations equal to thecenter value are included in both subsets.

Usage

dualSD(x, na.rm = FALSE, nmin = 10, center = xbar)

Arguments

Details

The purpose of dual SDs is to describe variability for asymmetric distributions. Symmetric distributions are also handled, though slightly less efficiently than a single SD does.

Value

a 2-vector of SDs with namesbottom andtop

Author(s)

Frank Harrell

See Also

pMedian()

Examples

set.seed(1)x <- rnorm(20000)sd(x)dualSD(x)y <- exp(x)s1 <- sd(y)s2 <- dualSD(y)s1s2quantile(y, c(0.025, 0.975))mean(y) + 1.96 * c(-1, 1) * s1mean(y) + 1.96 * c(- s2['bottom'], s2['top'])c(mean=mean(y), pseudomedian=pMedian(y), median=median(y))

ebpcomp

Description

Computation of Coordinates of Extended Box Plots Elements

Usage

ebpcomp(x, qref = c(0.5, 0.25, 0.75), probs = c(0.05, 0.125, 0.25, 0.375))

Arguments

Details

For an extended box plots computes all the elements needed for plotting it. This is typically used when adding to aggplot2 plot.

Value

list with elementssegments,lines,points,points2

Author(s)

Frank Harrell

Examples

ebpcomp(1:1000)

ecdfSteps

Description

Compute Coordinates of an Empirical Distribution Function

Usage

ecdfSteps(x, extend)

Arguments

Details

For a numeric vector uses the R built-inecdf function to computecoordinates of the ECDF, with extension slightly below and above therange ofx by default. This is useful forggplot2 where the ECDF may need to be transformed. The returned object is suitable for creating stratified statistics usingdata.table and other methods.

Value

a list with componentsx andy

Author(s)

Frank Harrell

See Also

stats::ecdf()

Examples

ecdfSteps(0:10)## Not run: # Use data.table for obtaining ECDFs by country and regionw <- d[, ecdfSteps(z, extend=c(1,11)), by=.(country, region)]  # d is a DT# Use ggplot2 to make one graph with multiple regions' ECDFs# and use faceting for countriesggplot(w, aes(x, y, color=region)) + geom_step() +       facet_wrap(~ country)## End(Not run)

Multicolumn Formating

Description

Expands the width either supercolumns or the subcolumns so that thethe sum of the supercolumn widths is the same as the sum of thesubcolumn widths.

Usage

equalBins(widths, subwidths)

Arguments

Details

This determins the correct subwidths of each of various columns in a tablefor printing. The correct width of the multicolumns is deterimed bysumming the widths of it subcolumns.

Value

widths of the the columns for a table.

Author(s)

Charles Dupont

See Also

nchar,stringDims

Examples

mcols <- c("Group 1", "Group 2")mwidth <- nchar(mcols, type="width")spancols <- c(3,3)ccols <- c("a", "deer", "ad", "cat", "help", "bob")cwidth <- nchar(ccols, type="width")subwidths <- partition.vector(cwidth, spancols)equalBins(mwidth, subwidths)

Plot Error Bars

Description

Add vertical error bars to an existing plot or makes a newplot with error bars.

Usage

errbar(x, y, yplus, yminus, cap=0.015, main = NULL,       sub=NULL, xlab=as.character(substitute(x)),       ylab=if(is.factor(x) || is.character(x)) ""           else as.character(substitute(y)),       add=FALSE, lty=1, type='p', ylim=NULL,       lwd=1, pch=16, errbar.col, Type=rep(1, length(y)),        ...)

Arguments

Details

errbar adds vertical error bars to an existing plot or makes a newplot with error bars. It can also make a horizontal error bar plotthat shows error bars for group differences as well as bars forgroups. For the latter type of plot, the lower x-axis scalecorresponds to group estimates and the upper scale corresponds todifferences. The spacings of the two scales are identical but thescale for differences has its origin shifted so that zero may beincluded. If at least one of the confidence intervals includes zero,a vertical dotted reference line at zero is drawn.

Author(s)

Charles Geyer, University of Chicago. Modified by Frank Harrell,Vanderbilt University, to handle missing data, to add the parametersadd andlty, and to implement horizontal charts with differences.

Examples

set.seed(1)x <- 1:10y <- x + rnorm(10)delta <- runif(10)errbar( x, y, y + delta, y - delta )# Show bootstrap nonparametric CLs for 3 group means and for# pairwise differences on same graphgroup <- sample(c('a','b','d'), 200, TRUE)y     <- runif(200) + .25*(group=='b') + .5*(group=='d')cla <- smean.cl.boot(y[group=='a'],B=100,reps=TRUE)  # usually B=1000a   <- attr(cla,'reps')clb <- smean.cl.boot(y[group=='b'],B=100,reps=TRUE)b   <- attr(clb,'reps')cld <- smean.cl.boot(y[group=='d'],B=100,reps=TRUE)d   <- attr(cld,'reps')a.b <- quantile(a-b,c(.025,.975))a.d <- quantile(a-d,c(.025,.975))b.d <- quantile(b-d,c(.025,.975))errbar(c('a','b','d','a - b','a - d','b - d'),       c(cla[1],clb[1],cld[1],cla[1]-clb[1],cla[1]-cld[1],clb[1]-cld[1]),       c(cla[3],clb[3],cld[3],a.b[2],a.d[2],b.d[2]),       c(cla[2],clb[2],cld[2],a.b[1],a.d[1],b.d[1]),       Type=c(1,1,1,2,2,2), xlab='', ylab='')

Escapes any characters that would have special meaning in a reqular expression.

Description

Escapes any characters that would have special meaning in a reqular expression.

Usage

escapeRegex(string)escapeBS(string)

Arguments

Details

escapeRegex will escape any characters that would havespecial meaning in a reqular expression. For any stringgrep(regexpEscape(string), string) will always be true.

escapeBS will escape any backslash ‘⁠\⁠’ in a string.

Value

The value of the string with any characters that would havespecial meaning in a reqular expression escaped.

Author(s)

Charles Dupont
Department of Biostatistics
Vanderbilt University

See Also

grep

Examples

string <- "this\\(system) {is} [full]."escapeRegex(string)escapeBS(string)

estSeqMarkovOrd

Description

Simulate Comparisons For Use in Sequential Markov Longitudinal Clinical Trial Simulations

Usage

estSeqMarkovOrd(  y,  times,  initial,  absorb = NULL,  intercepts,  parameter,  looks,  g,  formula,  ppo = NULL,  yprevfactor = TRUE,  groupContrast = NULL,  cscov = FALSE,  timecriterion = NULL,  coxzph = FALSE,  sstat = NULL,  rdsample = NULL,  maxest = NULL,  maxvest = NULL,  nsim = 1,  progress = FALSE,  pfile = "")

Arguments

Details

Simulates sequential clinical trials of longitudinal ordinal outcomes using a first-order Markov model. Looks are done sequentially after subject ID numbers given in the vectorlooks with the earliest possible look being after subject 2. At each look, a subject's repeated records are either all used or all ignored depending on the sequent ID number. For each true effect parameter value, simulation, and at each look, runs a function to compute the estimate of the parameter of interest along with its variance. For each simulation, data are first simulated for the last look, and these data are sequentially revealed for earlier looks. The user provides a functiong that has extra arguments specifying the true effect ofparameter the treatmentgroup expecting treatments to be coded 1 and 2.parameter is usually on the scale of a regression coefficient, e.g., a log odds ratio. Fitting is done using therms::lrm() function, unless non-proportional odds is allowed in which caseVGAM::vglm() is used. Iftimecriterion is specified, the function also, for the last data look only, computes the first time at which the criterion is satisfied for the subject or use the event time and event/censoring indicator computed bytimecriterion. The Cox/logrank chi-square statistic for comparing groups on the derived time variable is saved. Ifcoxzph=TRUE, thesurvival package correlation coefficientrho from the scaled partial residuals is also saved so that the user can later determine to what extent the Markov model resulted in the proportional hazards assumption being violated when analyzing on the time scale.vglm is accelerated by saving the first successful fit for the largest sample size and using its coefficients as starting value for furthervglm fits for any sample size for the same setting ofparameter.

Value

a data frame with number of rows equal to the product ofnsim, the length oflooks, and the length ofparameter, with variablessim,parameter,look,est (log odds ratio for group), andvest (the variance of the latter). Iftimecriterion is specified the data frame also containsloghr (Cox log hazard ratio for group),lrchisq (chi-square from Cox test for group), and ifcoxph=TRUE,phchisq, the chi-square for testing proportional hazards. The attributeetimefreq is also present iftimecriterion is present, and it probvides the frequency distribution of derived event times by group and censoring/event indicator. Ifsstat is given, the attributesstat is also present, and it contains an array with dimensions corresponding to simulations, parameter values within simulations,id, and a two-column subarray with columnsgroup andy, the latter being the summary measure computed by thesstat function. The returned data frame also has attributelrmcoef which are the last-look logistic regression coefficient estimates over thensim simulations and the parameter settings, and an attributefailures which is a data frame containing the variablesreason andfrequency cataloging the reasons for unsuccessful model fits.

Author(s)

Frank Harrell

See Also

gbayesSeqSim(),simMarkovOrd(),https://hbiostat.org/R/Hmisc/markov/

estSeqSim

Description

Simulate Comparisons For Use in Sequential Clinical Trial Simulations

Usage

estSeqSim(parameter, looks, gendat, fitter, nsim = 1, progress = FALSE)

Arguments

Details

Simulates sequential clinical trials. Looks are done sequentially at observation numbers given in the vectorlooks with the earliest possible look being at observation 2. For each true effect parameter value, simulation, and at each look, runs a function to compute the estimate of the parameter of interest along with its variance. For each simulation, data are first simulated for the last look, and these data are sequentially revealed for earlier looks. The user provides a functiongendat that given a true effect ofparameter and the two sample sizes (for treatment groups 1 and 2) returns a list with vectorsy1 andy2 containing simulated data. The user also provides a functionfitter with argumentsx (group indicator 0/1) andy (response variable) that returns a 2-vector containing the effect estimate and its variance.parameter is usually on the scale of a regression coefficient, e.g., a log odds ratio.

Value

a data frame with number of rows equal to the product ofnsim, the length oflooks, and the length ofparameter.

Author(s)

Frank Harrell

See Also

gbayesSeqSim(),simMarkovOrd(),estSeqMarkovOrd()

Examples

if (requireNamespace("rms", quietly = TRUE)) {  # Run 100 simulations, 5 looks, 2 true parameter values  # Total simulation time: 2s  lfit <- function(x, y) {  f <- rms::lrm.fit(x, y)    k <- length(coef(f))    c(coef(f)[k], vcov(f)[k, k])  }  gdat <- function(beta, n1, n2) {    # Cell probabilities for a 7-category ordinal outcome for the control group    p <- c(2, 1, 2, 7, 8, 38, 42) / 100    # Compute cell probabilities for the treated group    p2 <- pomodm(p=p, odds.ratio=exp(beta))    y1 <- sample(1 : 7, n1, p,  replace=TRUE)    y2 <- sample(1 : 7, n2, p2, replace=TRUE)    list(y1=y1, y2=y2)  }  set.seed(1)  est <- estSeqSim(c(0, log(0.7)), looks=c(50, 75, 95, 100, 200),                    gendat=gdat,                    fitter=lfit, nsim=100)  head(est)}

Flexible Event Chart for Time-to-Event Data

Description

Creates an event chart on the current graphics device. Also, allows userto plot legend on plot area or on separate page.Contains features useful for plotting data with time-to-event outcomesWhich arise in a variety of studiesincluding randomized clinical trials and non-randomized cohort studies.This function can use as input a matrix or a data frame, although greaterutility and ease of use will be seen with a data frame.

Usage

event.chart(data, subset.r = 1:dim(data)[1], subset.c = 1:dim(data)[2],           sort.by = NA, sort.ascending = TRUE,           sort.na.last = TRUE, sort.after.subset = TRUE,           y.var = NA, y.var.type = "n",           y.jitter = FALSE, y.jitter.factor = 1,           y.renum = FALSE, NA.rm = FALSE, x.reference = NA,           now = max(data[, subset.c], na.rm = TRUE),           now.line = FALSE, now.line.lty = 2,           now.line.lwd = 1, now.line.col = 1, pty = "m",           date.orig = c(1, 1, 1960), titl = "Event Chart",           y.idlabels = NA, y.axis = "auto",           y.axis.custom.at = NA, y.axis.custom.labels = NA,           y.julian = FALSE, y.lim.extend = c(0, 0),           y.lab = ifelse(is.na(y.idlabels), "", as.character(y.idlabels)),           x.axis.all = TRUE, x.axis = "auto",           x.axis.custom.at = NA, x.axis.custom.labels = NA,           x.julian = FALSE, x.lim.extend = c(0, 0), x.scale = 1,           x.lab = ifelse(x.julian, "Follow-up Time", "Study Date"),           line.by = NA, line.lty = 1, line.lwd = 1, line.col = 1,           line.add = NA, line.add.lty = NA,           line.add.lwd = NA, line.add.col = NA,           point.pch = 1:length(subset.c),           point.cex = rep(0.6, length(subset.c)),           point.col = rep(1, length(subset.c)),           point.cex.mult = 1., point.cex.mult.var = NA,           extra.points.no.mult = rep(NA, length(subset.c)),           legend.plot = FALSE, legend.location = "o", legend.titl = titl,           legend.titl.cex = 3, legend.titl.line = 1,           legend.point.at = list(x = c(5, 95), y = c(95, 30)),           legend.point.pch = point.pch,           legend.point.text = ifelse(rep(is.data.frame(data), length(subset.c)),                                      names(data[, subset.c]),                                      subset.c),           legend.cex = 2.5, legend.bty = "n",           legend.line.at = list(x = c(5, 95), y = c(20, 5)),           legend.line.text = names(table(as.character(data[, line.by]),                                          exclude = c("", "NA"))),           legend.line.lwd = line.lwd, legend.loc.num = 1,           ...)

Arguments