R : Copyright 2005, The R Foundation for Statistical Computing Version 2.1.1 (2005-06-20), ISBN 3-900051-07-0 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for a HTML browser interface to help. Type 'q()' to quit R. > ### *
> ### > attach(NULL, name = "CheckExEnv") > assign(".CheckExEnv", as.environment(2), pos = length(search())) # base > ## add some hooks to label plot pages for base and grid graphics > setHook("plot.new", ".newplot.hook") > setHook("persp", ".newplot.hook") > setHook("grid.newpage", ".gridplot.hook") > > assign("cleanEx", + function(env = .GlobalEnv) { + rm(list = ls(envir = env, all.names = TRUE), envir = env) + RNGkind("default", "default") + set.seed(1) + options(warn = 1) + delayedAssign("T", stop("T used instead of TRUE"), + assign.env = .CheckExEnv) + delayedAssign("F", stop("F used instead of FALSE"), + assign.env = .CheckExEnv) + sch <- search() + newitems <- sch[! sch %in% .oldSearch] + for(item in rev(newitems)) + eval(substitute(detach(item), list(item=item))) + missitems <- .oldSearch[! .oldSearch %in% sch] + if(length(missitems)) + warning("items ", paste(missitems, collapse=", "), + " have been removed from the search path") + }, + env = .CheckExEnv) > assign("..nameEx", "__{must remake R-ex/*.R}__", env = .CheckExEnv) # for now > assign("ptime", proc.time(), env = .CheckExEnv) > grDevices::postscript("Epi-Examples.ps") > assign("par.postscript", graphics::par(no.readonly = TRUE), env = .CheckExEnv) > options(contrasts = c(unordered = "contr.treatment", ordered = "contr.poly")) > options(warn = 1) > library('Epi') > > assign(".oldSearch", search(), env = .CheckExEnv) > assign(".oldNS", loadedNamespaces(), env = .CheckExEnv) > cleanEx(); ..nameEx <- "Lexis" > > ### * Lexis > > flush(stderr()); flush(stdout()) > > ### Name: Lexis > ### Title: Split follow-up time in cohort studies. > ### Aliases: Lexis > ### Keywords: manip > > ### ** Examples > > # A small bogus cohort > # > xcoh <- structure( list( id = c("A", "B", "C"), + birth = c("14/07/1952", "01/04/1954", "10/06/1987"), + entry = c("04/08/1965", "08/09/1972", "23/12/1991"), + exit = c("27/06/1997", "23/05/1995", "24/07/1998"), + fail = c(1, 0, 1) ), + .Names = c("id", "birth", "entry", "exit", "fail"), + row.names = c("1", "2", "3"), + class = "data.frame" ) > > # Convert the character dates into numerical variables (fractional years) > # > xcoh$bt <- cal.yr( xcoh$birth, format="%d/%m/%Y" ) > xcoh$en <- cal.yr( xcoh$entry, format="%d/%m/%Y" ) > xcoh$ex <- cal.yr( xcoh$exit , format="%d/%m/%Y" ) > > # See how it looks > # > xcoh id birth entry exit fail bt en ex 1 A 14/07/1952 04/08/1965 27/06/1997 1 1952.533 1965.589 1997.485 2 B 01/04/1954 08/09/1972 23/05/1995 0 1954.246 1972.686 1995.388 3 C 10/06/1987 23/12/1991 24/07/1998 1 1987.437 1991.974 1998.559 > > # Split time along one time-axis > # > Lexis( entry = en, + exit = ex, + fail = fail, + scale = 1, + origin = bt, + breaks = seq( 5, 40, 5 ), + include = list( bt, en, ex, id ), + data = xcoh ) Expand Entry Exit Fail Time bt en ex id 1 1 1965.589 1967.533 0 10 1952.533 1965.589 1997.485 A 2 1 1967.533 1972.533 0 15 1952.533 1965.589 1997.485 A 3 1 1972.533 1977.533 0 20 1952.533 1965.589 1997.485 A 4 1 1977.533 1982.533 0 25 1952.533 1965.589 1997.485 A 5 1 1982.533 1987.533 0 30 1952.533 1965.589 1997.485 A 6 1 1987.533 1992.533 0 35 1952.533 1965.589 1997.485 A 7 2 1972.686 1974.246 0 15 1954.246 1972.686 1995.388 B 8 2 1974.246 1979.246 0 20 1954.246 1972.686 1995.388 B 9 2 1979.246 1984.246 0 25 1954.246 1972.686 1995.388 B 10 2 1984.246 1989.246 0 30 1954.246 1972.686 1995.388 B 11 2 1989.246 1994.246 0 35 1954.246 1972.686 1995.388 B 12 3 1992.437 1997.437 0 5 1987.437 1991.974 1998.559 C 13 3 1997.437 1998.559 1 10 1987.437 1991.974 1998.559 C > > # Split time along two time-axes > # > ( x2 <- + Lexis( entry = en, + exit = ex, + fail = fail, + scale = 1, + origin = list( per=0, age=bt ), + breaks = list( per=seq(1900,2000,10), age=seq(0,80,5) ), + include = list( bt, en, ex, id ), + data = xcoh ) ) Expand Entry Exit Fail per age bt en ex id 1 1 1965.589 1967.533 0 1960 10 1952.533 1965.589 1997.485 A 2 1 1967.533 1970.000 0 1960 15 1952.533 1965.589 1997.485 A 3 1 1970.000 1972.533 0 1970 15 1952.533 1965.589 1997.485 A 4 1 1972.533 1977.533 0 1970 20 1952.533 1965.589 1997.485 A 5 1 1977.533 1980.000 0 1970 25 1952.533 1965.589 1997.485 A 6 1 1980.000 1982.533 0 1980 25 1952.533 1965.589 1997.485 A 7 1 1982.533 1987.533 0 1980 30 1952.533 1965.589 1997.485 A 8 1 1987.533 1990.000 0 1980 35 1952.533 1965.589 1997.485 A 9 1 1990.000 1992.533 0 1990 35 1952.533 1965.589 1997.485 A 10 1 1992.533 1997.485 1 1990 40 1952.533 1965.589 1997.485 A 11 2 1972.686 1974.246 0 1970 15 1954.246 1972.686 1995.388 B 12 2 1974.246 1979.246 0 1970 20 1954.246 1972.686 1995.388 B 13 2 1979.246 1980.000 0 1970 25 1954.246 1972.686 1995.388 B 14 2 1980.000 1984.246 0 1980 25 1954.246 1972.686 1995.388 B 15 2 1984.246 1989.246 0 1980 30 1954.246 1972.686 1995.388 B 16 2 1989.246 1990.000 0 1980 35 1954.246 1972.686 1995.388 B 17 2 1990.000 1994.246 0 1990 35 1954.246 1972.686 1995.388 B 18 2 1994.246 1995.388 0 1990 40 1954.246 1972.686 1995.388 B 19 3 1991.974 1992.437 0 1990 0 1987.437 1991.974 1998.559 C 20 3 1992.437 1997.437 0 1990 5 1987.437 1991.974 1998.559 C 21 3 1997.437 1998.559 1 1990 10 1987.437 1991.974 1998.559 C > > # Tabulate the cases and the person-years > # > tapply( x2$Fail, list( x2$age, x2$per ), sum ) 1960 1970 1980 1990 0 NA NA NA 0 5 NA NA NA 0 10 0 NA NA 1 15 0 0 NA NA 20 NA 0 NA NA 25 NA 0 0 NA 30 NA NA 0 NA 35 NA NA 0 0 40 NA NA NA 1 > tapply( x2$Exit - x2$Entry, list( x2$age, x2$per ), sum ) 1960 1970 1980 1990 0 NA NA NA 0.4633818 5 NA NA NA 5.0000000 10 1.943190 NA NA 1.1211493 15 2.467488 4.093088 NA NA 20 NA 10.000000 NA NA 25 NA 3.221081 6.778919 NA 30 NA NA 10.000000 NA 35 NA NA 3.221081 6.7789190 40 NA NA NA 6.0944547 > > > > cleanEx(); ..nameEx <- "Lexis.diagram" > > ### * Lexis.diagram > > flush(stderr()); flush(stdout()) > > ### Name: Lexis.diagram > ### Title: Plot a Lexis diagram > ### Aliases: Lexis.diagram > ### Keywords: hplot dplot > > ### ** Examples > > Lexis.diagram( entry.age = c(3,30,45), + risk.time = c(25,5,14), + birth.date = c(1970,1931,1925.7), + fail = c(TRUE,TRUE,FALSE) ) > LL <- Lexis.diagram( entry.age = sample( 0:50, 17, replace=TRUE ), + risk.time = sample( 5:40, 17, r=TRUE), + birth.date = sample( 1910:1980, 17, r=TRUE ), + fail = sample( 0:1, 17, r=TRUE ), + cex.fail = 1.1, + lwd.life = 2 ) > # Identify the persons' entry and exits > text( LL$exit.date, LL$exit.age, paste(1:nrow(LL)), col="red", font=2, adj=c(0,1) ) > text( LL$entry.date, LL$entry.age, paste(1:nrow(LL)), col="blue", font=2, adj=c(1,0) ) > data( nickel ) > attach( nickel ) > LL <- Lexis.diagram( age=c(10,100), date=c(1900,1990), + entry.age=age1st, exit.age=ageout, birth.date=dob, + fail=(icd %in% c(162,163)), lwd.life=1, + cex.fail=0.8, col.fail=c("green","red") ) > abline( v=1934, col="blue" ) > nickel[1:10,] id icd exposure dob age1st agein ageout 1 3 0 5.0 1889.019 17.4808 45.2273 92.9808 2 4 162 5.0 1885.978 23.1864 48.2684 63.2712 3 6 163 10.0 1881.255 25.2452 52.9917 54.1644 4 8 527 9.0 1886.340 24.7206 47.9067 69.6794 5 9 150 0.0 1879.500 29.9575 54.7465 76.8442 6 10 163 2.0 1889.915 21.2877 44.3314 62.5413 7 15 334 0.0 1890.500 23.2836 43.7465 62.0000 8 16 160 0.5 1874.332 50.3569 59.9149 65.5835 9 17 420 0.0 1909.500 15.4618 34.7465 50.5137 10 18 12 0.0 1892.500 24.1394 51.7465 57.5931 > LL[1:10,] entry.date entry.age exit.date exit.age birth.date risk.time fail 1 1906.500 17.4808 1982.000 92.9808 1889.019 75.5000 FALSE 2 1909.164 23.1864 1949.249 63.2712 1885.978 40.0848 TRUE 3 1906.500 25.2452 1935.419 54.1644 1881.255 28.9192 TRUE 4 1911.060 24.7206 1956.019 69.6794 1886.340 44.9588 FALSE 5 1909.457 29.9575 1956.344 76.8442 1879.500 46.8867 FALSE 6 1911.203 21.2877 1952.456 62.5413 1889.915 41.2536 TRUE 7 1913.784 23.2836 1952.500 62.0000 1890.500 38.7164 FALSE 8 1924.688 50.3569 1939.915 65.5835 1874.332 15.2266 FALSE 9 1924.962 15.4618 1960.014 50.5137 1909.500 35.0519 FALSE 10 1916.639 24.1394 1950.093 57.5931 1892.500 33.4537 FALSE > > > > cleanEx(); ..nameEx <- "Lexis.lines" > > ### * Lexis.lines > > flush(stderr()); flush(stdout()) > > ### Name: Lexis.lines > ### Title: Draw life lines in a Lexis diagram. > ### Aliases: Lexis.lines > ### Keywords: hplot dplot > > ### ** Examples > > Lexis.diagram( entry.age = c(3,30,45), + risk.time = c(25,5,14), + birth.date = c(1970,1931,1925.7), + fail = c(TRUE,TRUE,FALSE) ) > Lexis.lines( entry.age = sample( 0:50, 100, replace=TRUE ), + risk.time = sample( 5:40, 100, r=TRUE), + birth.date = sample( 1910:1980, 100, r=TRUE ), + fail = sample(0:1,100,r=TRUE), + cex.fail = 0.5, + lwd.life = 1 ) > > > > cleanEx(); ..nameEx <- "Life.lines" > > ### * Life.lines > > flush(stderr()); flush(stdout()) > > ### Name: Life.lines > ### Title: Compute dates/ages for life lines in a Lexis diagram > ### Aliases: Life.lines > ### Keywords: manip dplot > > ### ** Examples > > ( Life.lines( entry.age = c(3,30,45), + risk.time = c(25,5,14), + birth.date = c(1970,1931,1925.7) ) ) entry.date entry.age exit.date exit.age birth.date risk.time 1 1973.0 3 1998.0 28 1970.0 25 2 1961.0 30 1966.0 35 1931.0 5 3 1970.7 45 1984.7 59 1925.7 14 > > # Draw a Lexis diagram > Lexis.diagram() > > # Compute entry and exit age and date. > ( LL <- Life.lines( entry.age = c(3,30,45), + risk.time = c(25,5,14), + birth.date = c(1970,1931,1925.7) ) ) entry.date entry.age exit.date exit.age birth.date risk.time 1 1973.0 3 1998.0 28 1970.0 25 2 1961.0 30 1966.0 35 1931.0 5 3 1970.7 45 1984.7 59 1925.7 14 > segments( LL[,1], LL[,2], LL[,3], LL[,4] ) # Plot the life lines. > > # Compute entry and exit age and date, supplying a date variable > bd <- ( c(1970,1931,1925.7) - 1970 ) * 365.25 > class( bd ) <- "Date" > ( Life.lines( entry.age = c(3,30,45), + risk.time = c(25,5,14), + birth.date = bd ) ) entry.date entry.age exit.date exit.age birth.date risk.time 1 1972-12-31 1095.75 1998-01-01 10227.00 1970-01-01 9131.25 2 1961-01-05 10957.50 1966-01-05 12783.75 1931-01-06 1826.25 3 1970-09-18 16436.25 1984-09-17 21549.75 1925-09-18 5113.50 > > > > cleanEx(); ..nameEx <- "ROC" > > ### * ROC > > flush(stderr()); flush(stdout()) > > ### Name: ROC > ### Title: Function to compute and draw ROC-curves. > ### Aliases: ROC > ### Keywords: manip htest > > ### ** Examples > > x <- rnorm( 100 ) > z <- rnorm( 100 ) > w <- rnorm( 100 ) > tigol <- function( x ) 1 - ( 1 + exp( x ) )^(-1) > y <- rbinom( 100, 1, tigol( 0.3 + 3*x + 5*z + 7*w ) ) > ROC( form = y ~ x + z, plot="ROC" ) Warning in if (grid) abline(h = 0:10/10, v = 0:10/10, col = gray(0.9)) : the condition has length > 1 and only the first element will be used > > > > cleanEx(); ..nameEx <- "Relevel" > > ### * Relevel > > flush(stderr()); flush(stdout()) > > ### Name: Relevel > ### Title: Reorder and combine levels of a factor > ### Aliases: Relevel > ### Keywords: manip > > ### ** Examples > > ff <- factor( sample( letters[1:5], 100, replace=TRUE ) ) > table( ff, Relevel( ff, list( AB=1:2, "Dee"=4, c(3,5) ) ) ) ff AB Dee c+e a 13 0 0 b 25 0 0 c 0 0 19 d 0 26 0 e 0 0 17 > table( ff, rr=Relevel( ff, list( 5:4, Z=c("c","a") ), coll="-und-", first=FALSE ) ) rr ff b e-und-d Z a 0 0 13 b 25 0 0 c 0 0 19 d 0 26 0 e 0 17 0 > > > > cleanEx(); ..nameEx <- "S.typh" > > ### * S.typh > > flush(stderr()); flush(stdout()) > > ### Name: S.typh > ### Title: Salmonella Typhimurium outbreak 1996 in Denmark. > ### Aliases: S.typh > ### Keywords: datasets > > ### ** Examples > > data(S.typh) > > > > cleanEx(); ..nameEx <- "apc.frame" > > ### * apc.frame > > flush(stderr()); flush(stdout()) > > ### Name: apc.frame > ### Title: Produce an empty frame for display of parameter-estimates from > ### Age-Period-Cohort-models. > ### Aliases: apc.frame > ### Keywords: hplot > > ### ** Examples > > par( mar=c(4,4,1,4) ) > res <- + apc.frame( a.lab=seq(30,90,20), cp.lab=seq(1880,2000,30), r.lab=c(1,2,5,10,20,50), + a.tic=seq(30,90,10), cp.tic=seq(1880,2000,10), r.tic=c(1:10,1:5*10), + gap=27 ) > res cp.offset RR.fac 1763 5 > # What are the axes actually? > par(c("usr","xlog","ylog")) $usr [1] 30.0000000 237.0000000 -0.0679588 1.7669288 $xlog [1] FALSE $ylog [1] TRUE > # How to plot in the age-part: a point at (50,10) > points( 50, 10, pch=16, cex=2, col="blue" ) > # How to plot in the cohort-period-part: a point at (1960,0.3) > points( 1960-res[1], 0.3*res[2], pch=16, cex=2, col="red" ) > > > > graphics::par(get("par.postscript", env = .CheckExEnv)) > cleanEx(); ..nameEx <- "apc.lines" > > ### * apc.lines > > flush(stderr()); flush(stdout()) > > ### Name: apc.lines > ### Title: Plot APC-estimates in an APC-frame. > ### Aliases: apc.lines > ### Keywords: hplot > > ### ** Examples > > > > > cleanEx(); ..nameEx <- "bdendo" > > ### * bdendo > > flush(stderr()); flush(stdout()) > > ### Name: bdendo > ### Title: A case-control study of endometrial cancer > ### Aliases: bdendo > ### Keywords: datasets > > ### ** Examples > > data(bdendo) > > > > cleanEx(); ..nameEx <- "bdendo11" > > ### * bdendo11 > > flush(stderr()); flush(stdout()) > > ### Name: bdendo11 > ### Title: A 1:1 subset of the endometrial cancer case-control study > ### Aliases: bdendo11 > ### Keywords: datasets > > ### ** Examples > > data(bdendo11) > > > > cleanEx(); ..nameEx <- "births" > > ### * births > > flush(stderr()); flush(stdout()) > > ### Name: births > ### Title: Births in a London Hospital > ### Aliases: births > ### Keywords: datasets > > ### ** Examples > > data(births) > > > > cleanEx(); ..nameEx <- "blcaIT" > > ### * blcaIT > > flush(stderr()); flush(stdout()) > > ### Name: blcaIT > ### Title: Bladder cancer mortality in Italian males > ### Aliases: blcaIT > ### Keywords: datasets > > ### ** Examples > > data(blcaIT) > > > > cleanEx(); ..nameEx <- "brv" > > ### * brv > > flush(stderr()); flush(stdout()) > > ### Name: brv > ### Title: Bereavement in an elderly cohort > ### Aliases: brv > ### Keywords: datasets > > ### ** Examples > > data(brv) > > > > cleanEx(); ..nameEx <- "cal.yr" > > ### * cal.yr > > flush(stderr()); flush(stdout()) > > ### Name: cal.yr > ### Title: Function to convert character, factor and various date objects > ### into a number. > ### Aliases: cal.yr > ### Keywords: manip chron > > ### ** Examples > > birth <- c("14/07/1952","01/04/1954","10/06/1987","16/05/1990", + "01/01/1996","01/01/1997","01/01/1998","01/01/1999") > ( bt.yr <- cal.yr( birth, format="%d/%m/%Y" ) ) [1] 1952.533 1954.246 1987.437 1990.370 1995.999 1997.001 1998.000 1998.999 > # To get readable output: > data.frame( birth, bt.yr ) birth bt.yr 1 14/07/1952 1952.533 2 01/04/1954 1954.246 3 10/06/1987 1987.437 4 16/05/1990 1990.370 5 01/01/1996 1995.999 6 01/01/1997 1997.001 7 01/01/1998 1998.000 8 01/01/1999 1998.999 > > > > cleanEx(); ..nameEx <- "ccwc" > > ### * ccwc > > flush(stderr()); flush(stdout()) > > ### Name: ccwc > ### Title: Generate a nested case-control study > ### Aliases: ccwc > ### Keywords: datagen > > ### ** Examples > > # > # For the diet and heart dataset, create a nested case-control study > # using the age scale and matching on job > # > data(diet) > dietcc <- ccwc(doe, dox, chd, origin=dob, controls=2, data=diet, + include=energy, match=job) Sampling risk sets: ............................................. Warning in ccwc(doe, dox, chd, origin = dob, controls = 2, data = diet, : there were tied failure times > > > > cleanEx(); ..nameEx <- "ci.lin" > > ### * ci.lin > > flush(stderr()); flush(stdout()) > > ### Name: ci.lin > ### Title: Compute linear functions of parameters with s.e. > ### Aliases: ci.lin > ### Keywords: models regression > > ### ** Examples > > # Bogus data: > f <- factor( sample( letters[1:5], 200, replace=TRUE ) ) > g <- factor( sample( letters[1:3], 200, replace=TRUE ) ) > x <- rnorm( 200 ) > y <- 7 + as.integer( f ) * 3 + 2 * x + 1.7 * rnorm( 200 ) > > # Fit a simple model: > mm <- lm( y ~ x + f + g ) > ci.lin( mm ) Estimate StdErr z P 2.5% 97.5% (Intercept) 9.8699811 0.3766173 26.2069272 0.000000e+00 9.1318248 10.6081374 x 1.9675948 0.1322452 14.8783859 0.000000e+00 1.7083990 2.2267906 fb 3.0995925 0.4459483 6.9505649 3.638201e-12 2.2255500 3.9736351 fc 6.0475735 0.4508893 13.4125453 0.000000e+00 5.1638466 6.9313003 fd 9.4082193 0.4439493 21.1921035 0.000000e+00 8.5380947 10.2783439 fe 12.2280972 0.4633821 26.3887990 0.000000e+00 11.3198849 13.1363094 gb -0.2063660 0.3162740 -0.6524914 5.140842e-01 -0.8262516 0.4135195 gc -0.1516740 0.3317805 -0.4571517 6.475620e-01 -0.8019517 0.4986038 > ci.lin( mm, subset=3:6, diff=TRUE, fnam=FALSE ) Estimate StdErr z P 2.5% 97.5% b vs. c -2.947981 0.3953015 -7.45755 8.815171e-14 -3.722758 -2.173204 b vs. d -6.308627 0.3863870 -16.32723 0.000000e+00 -7.065931 -5.551322 b vs. e -9.128505 0.4193163 -21.76998 0.000000e+00 -9.950349 -8.306660 c vs. d -3.360646 0.3872787 -8.67759 0.000000e+00 -4.119698 -2.601594 c vs. e -6.180524 0.4252300 -14.53454 0.000000e+00 -7.013959 -5.347088 d vs. e -2.819878 0.4178048 -6.74927 1.485900e-11 -3.638760 -2.000995 > ci.lin( mm, subset=3:6, diff=TRUE, fnam=TRUE ) Estimate StdErr z P 2.5% 97.5% f b vs. c -2.947981 0.3953015 -7.45755 8.815171e-14 -3.722758 -2.173204 f b vs. d -6.308627 0.3863870 -16.32723 0.000000e+00 -7.065931 -5.551322 f b vs. e -9.128505 0.4193163 -21.76998 0.000000e+00 -9.950349 -8.306660 f c vs. d -3.360646 0.3872787 -8.67759 0.000000e+00 -4.119698 -2.601594 f c vs. e -6.180524 0.4252300 -14.53454 0.000000e+00 -7.013959 -5.347088 f d vs. e -2.819878 0.4178048 -6.74927 1.485900e-11 -3.638760 -2.000995 > ci.lin( mm, subset="f", diff=TRUE, fnam="f levels:" ) Estimate StdErr z P 2.5% f levels: a vs. b -3.099593 0.4459483 -6.950565 3.638201e-12 -3.973635 f levels: a vs. c -6.047573 0.4508893 -13.412545 0.000000e+00 -6.931300 f levels: a vs. d -9.408219 0.4439493 -21.192103 0.000000e+00 -10.278344 f levels: a vs. e -12.228097 0.4633821 -26.388799 0.000000e+00 -13.136309 f levels: b vs. c -2.947981 0.3953015 -7.457550 8.815171e-14 -3.722758 f levels: b vs. d -6.308627 0.3863870 -16.327225 0.000000e+00 -7.065931 f levels: b vs. e -9.128505 0.4193163 -21.769975 0.000000e+00 -9.950349 f levels: c vs. d -3.360646 0.3872787 -8.677590 0.000000e+00 -4.119698 f levels: c vs. e -6.180524 0.4252300 -14.534542 0.000000e+00 -7.013959 f levels: d vs. e -2.819878 0.4178048 -6.749270 1.485900e-11 -3.638760 97.5% f levels: a vs. b -2.225550 f levels: a vs. c -5.163847 f levels: a vs. d -8.538095 f levels: a vs. e -11.319885 f levels: b vs. c -2.173204 f levels: b vs. d -5.551322 f levels: b vs. e -8.306660 f levels: c vs. d -2.601594 f levels: c vs. e -5.347088 f levels: d vs. e -2.000995 > ci.lin( mm, subset="g", diff=TRUE, fnam="gee!:", vcov=TRUE ) > > # Use character defined subset to get ALL contrasts: > ci.lin( mm, subset="f", diff=TRUE ) Estimate StdErr z P 2.5% 97.5% a vs. b -3.099593 0.4459483 -6.950565 3.638201e-12 -3.973635 -2.225550 a vs. c -6.047573 0.4508893 -13.412545 0.000000e+00 -6.931300 -5.163847 a vs. d -9.408219 0.4439493 -21.192103 0.000000e+00 -10.278344 -8.538095 a vs. e -12.228097 0.4633821 -26.388799 0.000000e+00 -13.136309 -11.319885 b vs. c -2.947981 0.3953015 -7.457550 8.815171e-14 -3.722758 -2.173204 b vs. d -6.308627 0.3863870 -16.327225 0.000000e+00 -7.065931 -5.551322 b vs. e -9.128505 0.4193163 -21.769975 0.000000e+00 -9.950349 -8.306660 c vs. d -3.360646 0.3872787 -8.677590 0.000000e+00 -4.119698 -2.601594 c vs. e -6.180524 0.4252300 -14.534542 0.000000e+00 -7.013959 -5.347088 d vs. e -2.819878 0.4178048 -6.749270 1.485900e-11 -3.638760 -2.000995 > > > > cleanEx(); ..nameEx <- "ci.pd" > > ### * ci.pd > > flush(stderr()); flush(stdout()) > > ### Name: ci.pd > ### Title: Compute confidence limits for a difference of two independent > ### proportions. > ### Aliases: ci.pd > ### Keywords: distribution htest > > ### ** Examples > > ( a <- matrix( sample( 10:40, 4 ), 2, 2 ) ) [,1] [,2] [1,] 18 26 [2,] 21 35 > ci.pd( a ) 2 by 2 table analysis: ------------------------------------------------------ Outcome : Col 1 Comparing : Row 1 vs. Row 2 Col 1 Col 2 P(Col 1) 95% conf. interval Row 1 18 21 0.4615 0.3135 0.6167 Row 2 26 35 0.4262 0.3090 0.5524 95% conf. interval Relative Risk: 1.0828 0.6926 1.6929 Sample Odds Ratio: 1.1538 0.5140 2.5901 Conditional MLE Odds Ratio: 1.1522 0.4744 2.7940 Probability difference: 0.0353 -0.1567 0.2272 Exact P-value: 0.8369 Asymptotic P-value: 0.7287 ------------------------------------------------------ > twoby2( t(a) ) 2 by 2 table analysis: ------------------------------------------------------ Outcome : Col 1 Comparing : Row 1 vs. Row 2 Col 1 Col 2 P(Col 1) 95% conf. interval Row 1 18 21 0.4615 0.3135 0.6167 Row 2 26 35 0.4262 0.3090 0.5524 95% conf. interval Relative Risk: 1.0828 0.6926 1.6929 Sample Odds Ratio: 1.1538 0.5140 2.5901 Conditional MLE Odds Ratio: 1.1522 0.4744 2.7940 Probability difference: 0.0353 -0.1567 0.2272 Exact P-value: 0.8369 Asymptotic P-value: 0.7287 ------------------------------------------------------ > prop.test( t(a) ) 2-sample test for equality of proportions with continuity correction data: t(a) X-squared = 0.0197, df = 1, p-value = 0.8883 alternative hypothesis: two.sided 95 percent confidence interval: -0.1854080 0.2560260 sample estimates: prop 1 prop 2 0.4615385 0.4262295 > ( A <- array( sample( 10:40, 20 ), dim=c(2,2,5) ) ) , , 1 [,1] [,2] [1,] 16 37 [2,] 36 28 , , 2 [,1] [,2] [1,] 26 15 [2,] 11 14 , , 3 [,1] [,2] [1,] 25 39 [2,] 18 19 , , 4 [,1] [,2] [1,] 23 40 [2,] 27 22 , , 5 [,1] [,2] [1,] 24 31 [2,] 12 35 > print( ci.pd( A ) ) prob.diff lo hi 16/(16+36) - 37/(37+28) = -0.26153846 -0.5112809 0.002272471 26/(26+11) - 15/(15+14) = 0.18546132 -0.3469607 0.360161060 25/(25+18) - 39/(39+19) = -0.09101844 -0.5947699 0.136473082 23/(23+27) - 40/(40+22) = -0.18516129 -0.5833516 0.028132737 24/(24+12) - 31/(31+35) = 0.19696970 -0.2870718 0.323288012 > > > > cleanEx(); ..nameEx <- "diet" > > ### * diet > > flush(stderr()); flush(stdout()) > > ### Name: diet > ### Title: Diet and heart data > ### Aliases: diet > ### Keywords: datasets > > ### ** Examples > > data(diet) > # Illustrate the follow-up in a Lexis diagram > Lexis.diagram( age=c(30,75), date=c(1965,1990), + entry.date=cal.yr(doe), exit.date=cal.yr(dox), birth.date=cal.yr(dob), + fail=(fail>0), pch.fail=c(NA,16), col.fail=c(NA,"red"), cex.fail=1.0, + data=diet ) > > > > > cleanEx(); ..nameEx <- "ewrates" > > ### * ewrates > > flush(stderr()); flush(stdout()) > > ### Name: ewrates > ### Title: Rates of lung and nasal cancer mortality, and total mortality. > ### Aliases: ewrates > ### Keywords: datasets > > ### ** Examples > > data(ewrates) > str(ewrates) `data.frame': 150 obs. of 5 variables: $ year : num 1931 1931 1931 1931 1931 ... $ age : num 10 15 20 25 30 35 40 45 50 55 ... $ lung : num 1 2 6 14 30 68 149 274 431 586 ... $ nasal: num 0 0 0 0 1 1 3 5 10 15 ... $ other: num 1269 2201 3116 3024 3188 ... > > > > cleanEx(); ..nameEx <- "gmortDK" > > ### * gmortDK > > flush(stderr()); flush(stdout()) > > ### Name: gmortDK > ### Title: Population mortality rates for Denmark in 5-years age groups. > ### Aliases: gmortDK > ### Keywords: datasets > > ### ** Examples > > data(gmortDK) > > > > cleanEx(); ..nameEx <- "lep" > > ### * lep > > flush(stderr()); flush(stdout()) > > ### Name: lep > ### Title: An unmatched case-control study of leprosy incidence > ### Aliases: lep > ### Keywords: datasets > > ### ** Examples > > data(lep) > > > > cleanEx(); ..nameEx <- "lungDK" > > ### * lungDK > > flush(stderr()); flush(stdout()) > > ### Name: lungDK > ### Title: Male lung cancer incidence in Denmark > ### Aliases: lungDK > ### Keywords: datasets > > ### ** Examples > > data( lungDK ) > # Draw a Lexis diagram and show the number of cases in it. > attach( lungDK ) > Lexis.diagram( age=c(40,90), date=c(1943,1993), coh.grid=TRUE ) > text( Px, Ax, paste( D ), cex=0.7 ) > > > > cleanEx(); ..nameEx <- "mh" > > ### * mh > > flush(stderr()); flush(stdout()) > > ### Name: mh > ### Title: Mantel-Haenszel analyses of cohort and case-control studies > ### Aliases: mh > ### Keywords: htest > > ### ** Examples > > # If d and y are 3-way tables of cases and person-years > # observation formed by tabulation by two confounders > # (named "C1" and "C2") an exposure of interest ("E"), > # the following command will calculate an overall > # Mantel-Haenszel comparison of the first two exposure > # groups. > # > # Generate some bogus data > dnam <- list( E=c("low","medium","high"), C1=letters[1:2], C2=LETTERS[1:4] ) > d <- array( sample( 2:80, 24), + dimnames=dnam, dim=sapply( dnam, length ) ) > y <- array( abs( rnorm( 24, 227, 50 ) ), + dimnames=dnam, dim=sapply( dnam, length ) ) > mh(d, y, compare="E") $groups [1] "E" "low" "medium" $control [1] "C1" "C2" $type [1] "Rate ratio" $q [1] 161.2130 $r [1] 161.8798 $u [1] -0.6667823 $v [1] 158.5576 $ratio [1] 0.995881 $se.log.ratio [1] 0.0779466 $cl.lower [1] 0.8760446 $cl.upper [1] 1.13211 $chisq [1] 0.002804020 $p.value [1] 0.9577694 attr(,"class") [1] "mh" > # > # Or, if exposure levels named "low" and "high" are to be > # compared and these are not the first two levels of E : > # > mh(d, y, compare="E", levels=c("low", "high")) $groups [1] "E" NA NA $control [1] "C1" "C2" $type [1] "Rate ratio" $q [1] 176.1455 $r [1] 185.4986 $u [1] -9.353107 $v [1] 179.5288 $ratio [1] 0.9495786 $se.log.ratio [1] 0.07412436 $cl.lower [1] 0.840582 $cl.upper [1] 1.072708 $chisq [1] 0.4872788 $p.value [1] 0.4851437 attr(,"class") [1] "mh" > # > # If we wish to carry out an analysis which controls for C1, > # but examines the results at each level of C2: > # > mh(d, y, compare="E", by="C2") $groups [1] "E" "low" "medium" $control [1] "C1" $type [1] "Rate ratio" $q A B C D 43.80598 37.02829 41.11446 39.26424 $r A B C D 27.92845 33.60032 34.67521 65.67578 $u A B C D 15.877534 3.427969 6.439256 -26.411541 $v A B C D 34.38687 35.20816 38.27038 50.69214 $ratio A B C D 1.5685075 1.1020219 1.1857020 0.5978496 $se.log.ratio A B C D 0.1676510 0.1682221 0.1638417 0.1402068 $cl.lower A B C D 1.1904878 0.8356427 0.9055972 0.4747173 $cl.upper A B C D 2.066561 1.453315 1.552444 0.752920 $chisq A B C D 7.3311725 0.3337571 1.0834494 13.7609000 $p.value A B C D 0.0067768832 0.5634551239 0.2979271878 0.0002076126 attr(,"class") [1] "mh" > # > # It is also possible to look at rate ratios for every > # combination of C1 and C2 : > # > mh(d, y, compare="E", by=c("C1", "C2")) $groups [1] "E" "low" "medium" $control NULL $type [1] "Rate ratio" $q C2 C1 A B C D a 8.192843 34.241915 22.31777 10.73253 b 35.613139 2.786376 18.79669 28.53171 $r C2 C1 A B C D a 19.45554 25.030660 14.44625 42.80240 b 8.47291 8.569663 20.22895 22.87338 $u C2 C1 A B C D a -11.26270 9.211255 7.871517 -32.069865 b 27.14023 -5.783287 -1.432261 5.658324 $v C2 C1 A B C D a 12.38709 29.736043 18.65788 24.50017 b 21.99978 5.472112 19.61250 26.19197 $ratio C2 C1 A B C D a 0.421106 1.3679989 1.5448829 0.2507461 b 4.203177 0.3251442 0.9291975 1.2473759 $se.log.ratio C2 C1 A B C D a 0.2787698 0.1862629 0.2405626 0.2309401 b 0.2700149 0.4787136 0.2271115 0.2003342 $cl.lower C2 C1 A B C D a 0.2662269 1.0069981 1.0400386 0.1714991 b 2.6958287 0.1479472 0.6395445 0.8971984 $cl.upper C2 C1 A B C D a 0.6660868 1.8584156 2.294783 0.3666118 b 6.5533447 0.7145711 1.350036 1.7342282 $chisq C2 C1 A B C D a 10.24036 2.853346 3.3208905 41.978332 b 33.48180 6.112156 0.1045950 1.222383 $p.value C2 C1 A B C D a 1.374007e-03 0.09118397 0.06840477 9.229040e-11 b 7.193413e-09 0.01342552 0.74638359 2.688934e-01 attr(,"class") [1] "mh" > # > # If dimensions and levels of the table are unnamed, they must > # be referred to by number. > # > > > > cleanEx(); ..nameEx <- "mortDK" > > ### * mortDK > > flush(stderr()); flush(stdout()) > > ### Name: mortDK > ### Title: Population mortality rates for Denmark in 1-year age-classes. > ### Aliases: mortDK > ### Keywords: datasets > > ### ** Examples > > data(mortDK) > > > > cleanEx(); ..nameEx <- "ncut" > > ### * ncut > > flush(stderr()); flush(stdout()) > > ### Name: ncut > ### Title: Function to group a variable in intervals. > ### Aliases: ncut > ### Keywords: manip > > ### ** Examples > > br <- c(-2,0,1,2.5) > x <- c( rnorm( 10 ), br, -3, 3 ) > cbind( x, l=ncut( x, breaks=br, type="l" ), + m=ncut( x, breaks=br, type="m" ), + r=ncut( x, breaks=br, type="r" ) )[order(x),] x l m r [1,] -3.0000000 NA NA -2.0 [2,] -2.0000000 -2.0 -2.00 -2.0 [3,] -0.8356286 -2.0 -1.00 0.0 [4,] -0.8204684 -2.0 -1.00 0.0 [5,] -0.6264538 -2.0 -1.00 0.0 [6,] -0.3053884 -2.0 -1.00 0.0 [7,] 0.0000000 0.0 0.00 0.0 [8,] 0.1836433 0.0 0.50 1.0 [9,] 0.3295078 0.0 0.50 1.0 [10,] 0.4874291 0.0 0.50 1.0 [11,] 0.5757814 0.0 0.50 1.0 [12,] 0.7383247 0.0 0.50 1.0 [13,] 1.0000000 1.0 1.00 1.0 [14,] 1.5952808 1.0 1.75 2.5 [15,] 2.5000000 2.5 2.50 2.5 [16,] 3.0000000 2.5 NA NA > x <- rnorm( 200 ) > plot( x, ncut( x, breaks=br, type="l" ), pch=16, col="blue", ylim=range(x) ) > abline( 0, 1 ) > abline( v=br ) > points( x, ncut( x, breaks=br, type="r" ), pch=16, col="red" ) > points( x, ncut( x, breaks=br, type="m" ), pch=16, col="green" ) > > > > cleanEx(); ..nameEx <- "nice" > > ### * nice > > flush(stderr()); flush(stdout()) > > ### Name: nice > ### Title: Nice breakpoints > ### Aliases: nice > ### Keywords: manip > > ### ** Examples > > nice( exp( rnorm( 100 ) ), log=TRUE ) [1] 0.1 0.2 0.5 1.0 2.0 5.0 10.0 20.0 > > > > cleanEx(); ..nameEx <- "nickel" > > ### * nickel > > flush(stderr()); flush(stdout()) > > ### Name: nickel > ### Title: A Cohort of Nickel Smelters in South Wales > ### Aliases: nickel > ### Keywords: datasets > > ### ** Examples > > data(nickel) > str(nickel) `data.frame': 679 obs. of 7 variables: $ id : num 3 4 6 8 9 10 15 16 17 18 ... $ icd : num 0 162 163 527 150 163 334 160 420 12 ... $ exposure: num 5 5 10 9 0 2 0 0.5 0 0 ... $ dob : num 1889 1886 1881 1886 1880 ... $ age1st : num 17.5 23.2 25.2 24.7 30.0 ... $ agein : num 45.2 48.3 53.0 47.9 54.7 ... $ ageout : num 93.0 63.3 54.2 69.7 76.8 ... > > > > cleanEx(); ..nameEx <- "pctab" > > ### * pctab > > flush(stderr()); flush(stdout()) > > ### Name: pctab > ### Title: Create percentages in a table > ### Aliases: pctab > ### Keywords: manip methods array > > ### ** Examples > > Aye <- sample( c("Yes","Si","Oui"), 177, replace=TRUE ) > Bee <- sample( c("Hum","Buzz"), 177, replace=TRUE ) > Sea <- sample( c("White","Black","Red","Dead"), 177, replace=TRUE ) > A <- table( Aye, Bee, Sea ) > A , , Sea = Black Bee Aye Buzz Hum Oui 2 11 Si 8 13 Yes 7 6 , , Sea = Dead Bee Aye Buzz Hum Oui 5 10 Si 10 10 Yes 7 6 , , Sea = Red Bee Aye Buzz Hum Oui 5 6 Si 11 5 Yes 7 7 , , Sea = White Bee Aye Buzz Hum Oui 7 9 Si 4 12 Yes 3 6 > ftable( pctab( A ) ) , , Sea = Black Bee Aye Buzz Hum Oui 10.5 30.6 Si 24.2 32.5 Yes 29.2 24.0 , , Sea = Dead Bee Aye Buzz Hum Oui 26.3 27.8 Si 30.3 25.0 Yes 29.2 24.0 , , Sea = Red Bee Aye Buzz Hum Oui 26.3 16.7 Si 33.3 12.5 Yes 29.2 28.0 , , Sea = White Bee Aye Buzz Hum Oui 36.8 25 Si 12.1 30 Yes 12.5 24 , , Sea = All Bee Aye Buzz Hum Oui 100 100 Si 100 100 Yes 100 100 , , Sea = N Bee Aye Buzz Hum Oui 19 36 Si 33 40 Yes 24 25 Sea Black Dead Red White All N Aye Bee Oui Buzz 10.52632 26.31579 26.31579 36.84211 100.00000 19.00000 Hum 30.55556 27.77778 16.66667 25.00000 100.00000 36.00000 Si Buzz 24.24242 30.30303 33.33333 12.12121 100.00000 33.00000 Hum 32.50000 25.00000 12.50000 30.00000 100.00000 40.00000 Yes Buzz 29.16667 29.16667 29.16667 12.50000 100.00000 24.00000 Hum 24.00000 24.00000 28.00000 24.00000 100.00000 25.00000 > ftable( pctab( addmargins( A, 1 ), 3 ) ) , , Sea = Black Bee Aye Buzz Hum Oui 10.5 30.6 Si 24.2 32.5 Yes 29.2 24.0 Sum 22.4 29.7 , , Sea = Dead Bee Aye Buzz Hum Oui 26.3 27.8 Si 30.3 25.0 Yes 29.2 24.0 Sum 28.9 25.7 , , Sea = Red Bee Aye Buzz Hum Oui 26.3 16.7 Si 33.3 12.5 Yes 29.2 28.0 Sum 30.3 17.8 , , Sea = White Bee Aye Buzz Hum Oui 36.8 25.0 Si 12.1 30.0 Yes 12.5 24.0 Sum 18.4 26.7 , , Sea = All Bee Aye Buzz Hum Oui 100 100 Si 100 100 Yes 100 100 Sum 100 100 , , Sea = N Bee Aye Buzz Hum Oui 19 36 Si 33 40 Yes 24 25 Sum 76 101 Sea Black Dead Red White All N Aye Bee Oui Buzz 10.52632 26.31579 26.31579 36.84211 100.00000 19.00000 Hum 30.55556 27.77778 16.66667 25.00000 100.00000 36.00000 Si Buzz 24.24242 30.30303 33.33333 12.12121 100.00000 33.00000 Hum 32.50000 25.00000 12.50000 30.00000 100.00000 40.00000 Yes Buzz 29.16667 29.16667 29.16667 12.50000 100.00000 24.00000 Hum 24.00000 24.00000 28.00000 24.00000 100.00000 25.00000 Sum Buzz 22.36842 28.94737 30.26316 18.42105 100.00000 76.00000 Hum 29.70297 25.74257 17.82178 26.73267 100.00000 101.00000 > round( ftable( pctab( addmargins( A, 1 ), 3 ), row.vars=3 ), 1) , , Sea = Black Bee Aye Buzz Hum Oui 10.5 30.6 Si 24.2 32.5 Yes 29.2 24.0 Sum 22.4 29.7 , , Sea = Dead Bee Aye Buzz Hum Oui 26.3 27.8 Si 30.3 25.0 Yes 29.2 24.0 Sum 28.9 25.7 , , Sea = Red Bee Aye Buzz Hum Oui 26.3 16.7 Si 33.3 12.5 Yes 29.2 28.0 Sum 30.3 17.8 , , Sea = White Bee Aye Buzz Hum Oui 36.8 25.0 Si 12.1 30.0 Yes 12.5 24.0 Sum 18.4 26.7 , , Sea = All Bee Aye Buzz Hum Oui 100 100 Si 100 100 Yes 100 100 Sum 100 100 , , Sea = N Bee Aye Buzz Hum Oui 19 36 Si 33 40 Yes 24 25 Sum 76 101 Aye Oui Si Yes Sum Bee Buzz Hum Buzz Hum Buzz Hum Buzz Hum Sea Black 10.5 30.6 24.2 32.5 29.2 24.0 22.4 29.7 Dead 26.3 27.8 30.3 25.0 29.2 24.0 28.9 25.7 Red 26.3 16.7 33.3 12.5 29.2 28.0 30.3 17.8 White 36.8 25.0 12.1 30.0 12.5 24.0 18.4 26.7 All 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 N 19.0 36.0 33.0 40.0 24.0 25.0 76.0 101.0 > > > > cleanEx(); ..nameEx <- "plotEst" > > ### * plotEst > > flush(stderr()); flush(stdout()) > > ### Name: plotEst > ### Title: Plot estimates with confidence limits > ### Aliases: plotEst pointsEst linesEst > ### Keywords: hplot models > > ### ** Examples > > # Bogus data and a linear model > f <- factor( sample( letters[1:5], 100, replace=TRUE ) ) > x <- rnorm( 100 ) > y <- 5 + 2 * as.integer( f ) + 0.8 * x + rnorm(100) * 2 > m1 <- lm( y ~ f ) > > # Produce some confidence intervals for contrast to first level > ( cf <- summary( m1 )$coef[2:5,1:2] %*% rbind( c(1,1,1), 1.96*(c(0,-1,1) ) ) ) [,1] [,2] [,3] fb 1.751096 0.2119139 3.290279 fc 3.997887 2.3776946 5.618079 fd 5.426593 3.8975697 6.955615 fe 8.287793 6.6293347 9.946251 > > # Plots with increasing amount of bells and whistles > par( mfcol=c(3,2), mar=c(3,3,2,1) ) > plotEst( cf ) > plotEst( cf, grid=TRUE ) > plotEst( cf, grid=TRUE, cex=2, lwd=3 ) > plotEst( cf, grid=TRUE, cex=2, col.points="red", col.lines="green" ) > plotEst( cf, grid=TRUE, cex=2, col.points="red", col.lines="green", + xlog=TRUE, xtic=c(1:8), xlim=c(0.8,6) ) > rownames( cf )[1] <- "Contrast to fa:\n\n fb" > plotEst( cf, grid=TRUE, cex=2, col.points=rainbow(4), col.lines=rainbow(4), vref=1 ) > > > > graphics::par(get("par.postscript", env = .CheckExEnv)) > cleanEx(); ..nameEx <- "rateplot" > > ### * rateplot > > flush(stderr()); flush(stdout()) > > ### Name: rateplot > ### Title: Functions to plot rates from a table classified by age and > ### calendar time (period) > ### Aliases: rateplot Aplot Pplot Cplot > ### Keywords: hplot > > ### ** Examples > > data( blcaIT ) > attach(blcaIT) > > # Table of rates: > bl.rate <- tapply( D, list(age,period), sum ) / + tapply( Y, list(age,period), sum ) > bl.rate 1955 1960 1965 1970 1975 25 0.0000003000 3.000000e-07 1.000000e-07 0.0000004000 0.0000012000 30 0.0000017000 1.800000e-06 1.200000e-06 0.0000008000 0.0000009000 35 0.0000032000 3.100000e-06 3.500000e-06 0.0000042000 0.0000032000 40 0.0000104000 1.050000e-05 9.100000e-06 0.0000104000 0.0000127000 45 0.0000286000 2.520000e-05 2.610000e-05 0.0000304000 0.0000316000 50 0.0000664000 7.030000e-05 6.430000e-05 0.0000646000 0.0000847000 55 0.0001271000 1.339000e-04 1.459000e-04 0.0001464000 0.0001638000 60 0.0002011000 2.398000e-04 2.669000e-04 0.0002755000 0.0002853000 65 0.0002440000 3.316000e-04 4.211999e-04 0.0004777000 0.0005037000 70 0.0003281000 4.231000e-04 5.287000e-04 0.0006601000 0.0007464000 75 0.0004554001 4.793999e-04 6.204999e-04 0.0008464998 0.0010420998 > > # The four classical plots: > par( mfrow=c(2,2) ) > rateplot( bl.rate*10^6 ) > > # The labels on the vertical axis could be nicer: > rateplot( bl.rate*10^6, at=10^(-1:3), labels=c(0.1,1,10,100,1000) ) > > # More bells an whistles > par( mfrow=c(1,3), mar=c(3,3,1,1), oma=c(0,3,0,0), mgp=c(3,1,0)/1.6 ) > rateplot( bl.rate*10^6, ylab="", ann=TRUE, which=c("AC","PA","CA"), + at=10^(-1:3), labels=c(0.1,1,10,100,1000), + col=topo.colors(11), cex.ann=1.2 ) > > > > graphics::par(get("par.postscript", env = .CheckExEnv)) > cleanEx(); ..nameEx <- "stattable" > > ### * stattable > > flush(stderr()); flush(stdout()) > > ### Name: stat.table > ### Title: Tables of summary statistics > ### Aliases: stat.table print.stat.table > ### Keywords: iteration category > > ### ** Examples > > data(warpbreaks) > # A one-way table > stat.table(tension,list(count(),mean(breaks)),data=warpbreaks) ------------------------------- tension count() mean(breaks) ------------------------------- L 18 36.39 M 18 26.39 H 18 21.67 ------------------------------- > # The same table with informative labels > stat.table(index=list("Tension level"=tension),list(N=count(), + "mean number of breaks"=mean(breaks)),data=warpbreaks) -------------------------- Tension N mean level number of breaks -------------------------- L 18 36.39 M 18 26.39 H 18 21.67 -------------------------- > > # A two-way table > stat.table(index=list(tension,wool),mean(breaks),data=warpbreaks) -------------------------- ------wool------- tension A B -------------------------- L 44.56 28.22 M 24.00 28.78 H 24.56 18.78 -------------------------- > # The same table with margins over tension, but not wool > stat.table(index=list(tension,wool),mean(breaks),data=warpbreaks, + margins=c(TRUE, FALSE)) -------------------------- ------wool------- tension A B -------------------------- L 44.56 28.22 M 24.00 28.78 H 24.56 18.78 Total 31.04 25.26 -------------------------- > > # A table of column percentages > stat.table(list(tension,wool), percent(tension), data=warpbreaks) -------------------------- ------wool------- tension A B -------------------------- L 33.3 33.3 M 33.3 33.3 H 33.3 33.3 -------------------------- > # Cell percentages, with margins > stat.table(list(tension,wool),percent(tension,wool), margin=TRUE, + data=warpbreaks) ---------------------------------- ----------wool----------- tension A B Total ---------------------------------- L 16.7 16.7 33.3 M 16.7 16.7 33.3 H 16.7 16.7 33.3 Total 50.0 50.0 100.0 ---------------------------------- > > # A table with multiple statistics > # Note how each statistic has its own default precision > a <- stat.table(index=list(wool,tension), + contents=list(count(),mean(breaks),percent (wool)), + data=warpbreaks) > print(a) ------------------------------- ---------tension--------- wool L M H ------------------------------- A 9 9 9 44.56 24.00 24.56 50.0 50.0 50.0 B 9 9 9 28.22 28.78 18.78 50.0 50.0 50.0 ------------------------------- > # Print the percentages rounded to the nearest integer > print(a, digits=c(percent=0)) ------------------------------- ---------tension--------- wool L M H ------------------------------- A 9 9 9 44.56 24.00 24.56 50 50 50 B 9 9 9 28.22 28.78 18.78 50 50 50 ------------------------------- > > # An Epidemiological example with follow-up time > data(nickel) > str(nickel) `data.frame': 679 obs. of 7 variables: $ id : num 3 4 6 8 9 10 15 16 17 18 ... $ icd : num 0 162 163 527 150 163 334 160 420 12 ... $ exposure: num 5 5 10 9 0 2 0 0.5 0 0 ... $ dob : num 1889 1886 1881 1886 1880 ... $ age1st : num 17.5 23.2 25.2 24.7 30.0 ... $ agein : num 45.2 48.3 53.0 47.9 54.7 ... $ ageout : num 93.0 63.3 54.2 69.7 76.8 ... > > # Make a grouped version of the exposure variable > nickel$egr <- cut( nickel$exposure, breaks=c(0, 0.5, 5, 10, 100), right=FALSE ) > stat.table( egr, list( count(), percent(egr), mean( age1st ) ), data=nickel ) --------------------------------------------- egr count() percent(egr) mean(age1st) --------------------------------------------- [0,0.5) 290 42.7 25.87 [0.5,5) 235 34.6 27.29 [5,10) 111 16.3 27.67 [10,100) 43 6.3 27.19 --------------------------------------------- > > # Split the follow-up time by current age > nickel.ex <- + Lexis( entry=agein, exit=ageout, fail=icd %in% c(162,163), + origin=0, breaks=seq(0,100,20), + include=list( id, exposure, egr, age1st, icd ), data=nickel ) > str( nickel.ex ) `data.frame': 1471 obs. of 10 variables: $ Expand : num 1 1 1 2 2 3 4 4 5 5 ... $ Entry : num 45.2 60.0 80.0 48.3 60.0 ... $ Exit : num 60.0 80.0 93.0 60.0 63.3 ... $ Fail : num 0 0 0 0 1 1 0 0 0 0 ... $ Time : num 40 60 80 40 60 40 40 60 40 60 ... $ id : num 3 3 3 4 4 6 8 8 9 9 ... $ exposure: num 5 5 5 5 5 10 9 9 0 0 ... $ egr : Factor w/ 4 levels "[0,0.5)","[0.5,5)",..: 3 3 3 3 3 4 3 3 1 1 ... $ age1st : num 17.5 17.5 17.5 23.2 23.2 ... $ icd : num 0 0 0 162 162 163 527 527 150 150 ... > > # Table of rates > stat.table( Time, list( n=count(), N=count(id), D=sum(Fail), + "Rate/1000"=ratio(Fail,Exit-Entry,1000) ), + margin=1, data=nickel.ex ) ------------------------------------------ Time n N D Rate/1000 ------------------------------------------ 20 239 238 0.00 0.00 40 644 643 63.00 7.96 60 501 500 70.00 12.31 80 87 87 4.00 9.64 Total 1471 676 137.00 8.93 ------------------------------------------ > # Two-way table of rates and no. persons contributing > stat.table( list(age=Time, Exposure=egr), + list( N=count(id), D=sum(Fail), Y=sum((Exit-Entry)/1000), + Rate=ratio(Fail,Exit-Entry,1000) ), + margin=TRUE, data=nickel.ex ) ------------------------------------------------- -----------------Exposure----------------- age [0,0.5) [0.5,5) [5,10) [10,100) Total ------------------------------------------------- 20 142 83 13 0 238 0.00 0.00 0.00 NA 0.00 0.90 0.39 0.05 NA 1.33 0.00 0.00 0.00 NA 0.00 40 268 229 105 41 643 19.00 24.00 15.00 5.00 63.00 3.85 2.84 0.92 0.31 7.91 4.94 8.46 16.27 16.32 7.96 60 223 172 82 24 500 19.00 28.00 16.00 7.00 70.00 2.83 1.79 0.79 0.27 5.69 6.71 15.61 20.16 26.34 12.31 80 50 25 9 3 87 4.00 0.00 0.00 0.00 4.00 0.22 0.11 0.06 0.02 0.41 18.04 0.00 0.00 0.00 9.64 Total 288 235 111 43 676 42.00 52.00 31.00 12.00 137.00 7.80 5.13 1.83 0.59 15.35 5.38 10.14 16.99 20.22 8.93 ------------------------------------------------- > > > > cleanEx(); ..nameEx <- "tabplot" > > ### * tabplot > > flush(stderr()); flush(stdout()) > > ### Name: tabplot > ### Title: Graphical display of a 2-way contingency table > ### Aliases: tabplot > ### Keywords: hplot > > ### ** Examples > > b <- sample( letters[1:4], 300, replace=TRUE, prob=c(3,1,2,4)/10 ) > a <- rnorm( 300 ) - as.integer( factor( b ) ) / 8 > tb <- table( cut( a, -3:2 ), b ) > tabplot( tb ) > tabplot( tb, rowlabs="right", col=heat.colors ) > > # Very similar plots > ptb <- sweep( tb, 2, apply( tb, 2, sum ), "/" ) > par( mfrow=c(2,2) ) > barplot( ptb, space=0 ) > tabplot( tb, equal=TRUE, lwd=1 ) > tabplot( tb, equal=TRUE, lwd=1, rowlabs="l" ) > tabplot( tb, equal=FALSE, lwd=1, rowlabs="l" ) > > > > graphics::par(get("par.postscript", env = .CheckExEnv)) > cleanEx(); ..nameEx <- "thoro" > > ### * thoro > > flush(stderr()); flush(stdout()) > > ### Name: thoro > ### Title: Thorotrast Study > ### Aliases: thoro > ### Keywords: datasets > > ### ** Examples > > data(thoro) > str(thoro) `data.frame': 2470 obs. of 14 variables: $ id : num 1 2 3 4 5 6 7 8 9 10 ... $ sex : num 2 2 1 1 1 2 1 2 1 1 ... $ birthdat:Class 'Date' num [1:2470] -19501 -15398 -24553 -18862 -24497 ... $ contrast: num 1 1 1 1 1 1 1 1 1 1 ... $ injecdat:Class 'Date' num [1:2470] -11399 -9531 -12554 -12274 -11912 ... $ volume : num 22 80 10 10 10 20 10 40 34 10 ... $ exitdat :Class 'Date' num [1:2470] 2479 -1450 -3755 2669 -8990 ... $ exitstat: num 1 1 1 1 1 1 1 3 1 1 ... $ cause : num 2 8 2 2 14 14 3 NA 2 2 ... $ liverdat:Class 'Date' num [1:2470] -1450 -1450 NA 2669 NA ... $ liver : num 1 1 0 1 0 0 0 0 1 0 ... $ hepcc : num 0 0 0 0 0 0 0 0 0 0 ... $ chola : num 0 0 0 0 0 0 0 0 0 0 ... $ hmang : num 1 1 0 1 0 0 0 0 1 0 ... > > > > cleanEx(); ..nameEx <- "twoby2" > > ### * twoby2 > > flush(stderr()); flush(stdout()) > > ### Name: twoby2 > ### Title: Analysis of a two by two table > ### Aliases: twoby2 > ### Keywords: univar htest > > ### ** Examples > > Treat <- sample(c("A","B"), 50, rep=TRUE ) > Resp <- c("Yes","No")[1+rbinom(50,1,0.3+0.2*(Treat=="A"))] > twoby2( Treat, Resp ) 2 by 2 table analysis: ------------------------------------------------------ Outcome : No Comparing : A vs. B No Yes P(No) 95% conf. interval A 6 17 0.2609 0.1222 0.4723 B 10 17 0.3704 0.2122 0.5623 95% conf. interval Relative Risk: 0.7043 0.3024 1.6407 Sample Odds Ratio: 0.6000 0.1780 2.0223 Conditional MLE Odds Ratio: 0.6062 0.1454 2.3472 Probability difference: -0.1095 -0.3406 0.1466 Exact P-value: 0.5456 Asymptotic P-value: 0.4099 ------------------------------------------------------ > twoby2( table( Treat, Resp )[,2:1] ) # Comparison the other way round 2 by 2 table analysis: ------------------------------------------------------ Outcome : Yes Comparing : A vs. B Yes No P(Yes) 95% conf. interval A 17 6 0.7391 0.5277 0.8778 B 17 10 0.6296 0.4377 0.7878 95% conf. interval Relative Risk: 1.1739 0.8047 1.7126 Sample Odds Ratio: 1.6667 0.4945 5.6174 Conditional MLE Odds Ratio: 1.6497 0.4260 6.8765 Probability difference: 0.1095 -0.1466 0.3406 Exact P-value: 0.5456 Asymptotic P-value: 0.4099 ------------------------------------------------------ > > > > ### *