R Tips and Links

Helpful R Links

1. 2014 MSP Orientation
2. My R Class Notes    (Advanced)
3. Official R Site Search     RSeek (Google-type search for R related material)     Inside R (documentation for base R and packages)
4. Documentation (incl. Download)   (Hint: Try creating a bookmark to "C:\Program Files\R\rw2001\doc\html\rwin.html" in Windows; substitute your most recent version for "rw2001"; the link to "Search Engine and Keywords" is most helpful.)
5. Packages     Crantastic Package Page     R Example Graph Library
6. FAQ     R-help for asking questions   Bug reporting    
7. R Inferno: problems and solutions
8. R Studio: an integrated development environment for R
9. R color chart (pdf)
10. Intro to R
11. R Language Definition
12. Operator Precedence
13. R Data Import/Export
14. R Reference Guide (pdf, 450 pages, 12MB)
15. Evaluating the design of the R language (pdf)
16. Shalizi's guide to R Markdown
17. Simple R Markdown figure labels
18. cache_with_help.zip is an R markdown template that sets good cache settings, includes simple R markdown figure labels, and has markdown help at the bottom. To use it, unzip it in the same folder as "github_document". Then choose Template after File / New Markdown in RStudio.
19. aRrgh: a newcomer's (angry) guide to R (gripes from the CS community)
20. R inferno (traps and tips)
21. Advanced R Programming
22. Vectorization in R
23. R Data Import/Export
24. Big Data with R (and Python)
25. A good proposed Programming Style Guide
26. Areas of statistics: CRAN Task Views
27. Lumley's R Fundamentals (pdf)
28. R Reference Card
29. Peng's Debugging in R (pdf)
30. R Programming Resource Center
31. Mathematical annotation in plots ("plotmath")
32. R-News
33. R wiki
34. Spherula (R quick reference, scripts, book notes, etc.)
35. OmegaHat: interfaces to other languages
36. Package gdata has function read.xls() which can read Excel files. An alternative is the RODBC package.
37. Exchanging data between R and MS Windows apps (Excel, etc)
38. Rtools needed for building packages (choose "Download R for" and an operating system)
39. Nabble R forum
40. Kickstarting R
41. Using R for psychological research (personality-project)
42. York U R tips
43. Books
44. R Overview and Links
45. Theresa Scott's tutorial
46. O'Reilly free tutorial
47. R for Cats tutorial
48. Online courses in R from Statistics.com (expensive)
49. Missing data blog page
50. CUNY R Tutorial
51. Zoonek R Tutorial
52. Tips for Creating, Modifying, and Checking Data Frames
53. Practical Regression and Anova using R (pdf)
54. Non-Parametric Inference with R by Larry and Chad (pdf)
55. Example of running repeated measures in R to match SPSS (etc)
56. nlme() mixed models guide (pdf)

"Shiny": a package for creating interactive R applications with a web browser interface

1. Try it: To try shiny, just click on one of the links given here. Note that these apps run on a remote server. Normally you will develop, and perhaps deploy, your shiny apps from within your R session. Also note that I have turned on the optional "showcase" mode so that you can see the code below each app.
   • A shiny app that helps you get a feel for quantile normal plots   (hosted at shinyapps.io, RStudio's hosting service for Shiny apps. )
   • The same quantile normal app hosted at stat.cmu.edu   (ui.R)   (server.R)   [loads qqn.R and spssSkewKurtosis.R from http://www.stat.cmu.edu/hseltman/files/]
   • Many other statistics-related Shiny apps: Shimmer
2. Official Homepage: The Shiny website
3. Function Reference: Shiny functions
4. Official Tutorial: The shiny tutorial
5. Advanced Tips: Attali's Tips
6. Advanced Videos: Developer Conference
7. ShinyTex: shinyTex
8. Getting started:
   1. Start R in any directory
   2. Install shiny using install.packages("shiny") if it has not previously been installed.
   3. Run library("shiny") once per R session (or place this command in the .First() function).
   4. Place files named ui.R and server.R in the directory (or in another directory, e.g., named "foo"). You can create these files from scratch, or you may want to start with the files linked here, which implement a simple histogram vs. boxplot app.
   5. From R in your working directory, run runApp() (or runApp("foo") if the ui.R and server.R files were placed in another directory).
   6. Interact with the app in the supplied test window. In this window you can click "Open in Browser" to open the app in your default browser. Certain functions, such as downloading a file using downloadButton() work only in a browser, and not in the test window.
   7. Make whatever additions / changes you need to either file to change the app so that it does what you want it to do (see the official tutorial and Inside R for details).
   8. Note that you can change the ui.R and/or server.R files while the app is running, save the changes, and click "Reload" in your browser to see the effects of the changes without quitting the app. (Check the R console for possible error messages.)
   9. To quit the app, use the "escape" key in the R console window or close the test window or browser window.
9. Showcase mode: If you want your users to learn about how you made your app, you can turn on showcase mode. The simplest way to do this is to place a DESCRIPTION plain text file in the same directory as your ui.R and server.R files. Create your version by editing this example of a DESCRIPTION file. If you want additional text about the app to show along with your ui.R and server.R files, create a plain text file called Readme.md. Here is an example that demonstrates headers and spacing. Full formatting information for this file type is given at Wikipedia Markdown, but currently not all markdown features are implemented in Shiny.
10. Five ways to deploy your shiny app:
    1. Distribute your ui.R and server.R files to users. They only need to put them in a directory, load shiny in R (library("shiny"), and run an R command like runApp() or runApp("myShinyDirectory").
    2. Put your ui.R and server.R files in a github gist. Then users only need to load shiny and run an R command like runGist("myFirstShinyApp")
    3. Put your ui.R and server.R files in zip file (as a subdirectory) and put the zip file on your website. Then users only need to load shiny and run an R command like runUrl("myFirstShinyApp.zip")
    4. Let R Studio host your app on their server. You register with "shinyapps", load the "shinyapps" package in R, and run deployApp() to upload your app to the RStudio servers. After some processing, the function tells you the URL that corresponds to your app. Then your users only need to enter the URL in their browser.
    5. Set up your own server (advanced). Your users only need to enter the URL in their browser
11. In development: locator() equivalent
12. Shiny apps that do not isolate text inputs and then use an "action button" to initiate any response to the input text will initiate the response roughly after every other letter typed. A good compromise would be to change text inputs so that they "invalidate" (initiate a response) only when Enter or Tab is pressed. Based on information from ZJ in the Google Groups Shiny discussion forum, I suggest this solution. (The main disadvantage is that it must be manually re-done whenever the Shiny package is updated.)
    The core of the problem is that "shiny.js" sets the "key-up" event to trigger a response to text input. So my suggestion is to remove that bit of Javascript code from the shiny library file. Actually, the file "shiny.js" is a human readable, version of the code, but the real code is "minimized" in "shiny.min.js", and that is what must be fixed.
    In addition to just searching for the file, here is a way to find it. In R, run myPaths = .libPaths() followed by print(myPaths) to see the one or several places that packages live on your computer. Then run list.files(myPaths[1]), etc. until you find the "shiny" package. Within that shiny folder, you will find "shiny.min.js" in the shared subfolder of the "www" subfolder.
    For safety sake, first copy the file to a backup name. Then use a plain text editor to open the file and change

    subscribe:function(a,b){$(a).on("keyup.textInputBinding input.textInputBinding",function(a){b(!0)}),$(a).on("change.textInputBinding",function(a){b(!1)})}
    

    to

    subscribe:function(a,b){$(a).on("change.textInputBinding",function(a){b(!1)})}
    

    Save the file and the behavior of textInput will now be fixed.

Tips for Working in R

1. Use helpstart() to bring up help in a browser; the link to "Search Engine and Keywords" is most useful.
2. When looking at complex expressions, decode them by working from the inside out. E.g., here is a decomposition of some code to make a density plot of the product to two normal random variables using a sample of size 20. (Note the optional use of "tmp" to keep the lines so that they all use the same random numbers.)

        plot(density(apply(matrix(rnorm(40),20), 1, prod)))
        tmp = rnorm(40)   # a vector of 40 standard normal variates
        tmp
        matrix(tmp,20)  # put into a matrix of 20 rows and 40/20=2 columns
        apply(matrix(tmp,20), 1, prod)  # the 20 products
        density(apply(matrix(tmp,20), 1, prod)) # the density estimate
        plot(density(apply(matrix(tmp,20), 1, prod))) # the plot
       

3. Keep a text record of all working R commands needed to re-run your analysis. Ideally you should be able to source() the file and recreate your work, e.g. if your client finds an error in the data (which happens 98% of the time according to Seltman's Law of Data Analysis).
4. Under Linux, "ESS" (Emacs Speaks S) is usually the most efficient way to work. Briefly, you start emacs, then use "Alt-X R" to start R from within emacs. The ESS menu (with keyboard shortcuts) allows you to automatically run code that you write, among other features. The home page is ESS .
5. Write out TRUE and FALSE, because T and F can be redefined.
6. When defining and redefining columns of a data.frame, make liberal use of summary() and table(..., exclude=NULL) to verify that you accomplished what you tried to accomplish.
7. Important:Remember that table() ignores missing data; use table(..., exclude=NULL) to also see missing data.
8. You can use the .First function to automatically load libraries that you frequently use, or to perform other startup tasks. E.g.

    .First = function() {library(nlme); options(locatorBell=FALSE)} 

   will load the nlme library every time R starts up in the current directory. It also turns off the annoying sound associated with the locator() function.
9. Use this function to find large, unneeded objects that can be removed to free up space:

        sizes = function() {
          ob = objects(name=parent.frame())
          rslt = sapply(ob,function(x){object.size(get(x))})
          return(sort(rslt))
        }
        

10. Contrast testing in R (using C() or contrasts()), ignores your scaling, so although the t-values and p-values are correct, the estimates and standard errors (and any confidence intervals you construct from them) are incorrect. To do this correctly, use fit.contrast() in package "gmodels". E.g.

         x = factor(rep(LETTERS[1:3], each=20)); y = rnorm(60)
         m1 = aov(y~x)
         library(gmodels)
         cont = rbind(AvsBC = c(1, -1/2, -1/2), BvsC = c(0, 1, -1))
         fit.contrast(m1, "x", cont, conf.int=0.95)
         

11. If you are working on a public computer without write access to where most of R lives, you can still install packages to a private space (Windows example shown here, but it is similar on Linux). Make a directory you can write to, e.g., c:\\myPackages. In R, to install, e.g., package "mice" use

         install.packages("mice", "c:\\myPackages")
         

    Then each R session use

         library(mice, lib.loc="c:\\myPackages")
         

12. A system for documenting data analysis projects:

    Here is an idea for making R code that stores comments and results in a separate, readable file. This is especially nice when you might need to re-source() your code due to changes in the data or analysis (i.e., essentially always). Optionally, you can run reportLatex() after all of your report() commands to create a .tex file that is formatted better and incorporates graphical output (see below).

    (An alternative is sweave. Unlike sweave, report does not require you to understand latex, and it has only a single command to learn.)

    The code and more documentation are at report.R.

    Put these two lines near the top of your code:

         if (!exists("report")) source("http://www.stat.cmu.edu/~hseltman/files/report.R")
         report("Start of my report on project X", new=TRUE, prefix="myProjectX")
         

    This creates a file named "myProjectXYYYY-MM-DD.txt" with the quoted string in the first argument as the text at the top of the file. You can include "\n" in the first argument to write multiple lines in one call. You can optionally add the argument useTime=TRUE to include the time of creation along with the date in the file name if you want to keep multiple versions from the same day.

    Note that the variable "reportFileName" is created in your global environment and you should not delete this variable, at least while you are working on any one report.

    Now anytime in your code, you can include code of the form

     report(x) 

    or

     report(x, ..., z)

    to cause the value of x (or all of the variables x through z) to go to both the screen and the report file. This constructs the report on-the-fly as you work through your analysis. (If multiple arguments are used with report() and they are all strings or single numbers, then they are pasted together without any spaces between them (i.e., using sep="")).

    Note that you can manually erase errors from the report file using a text editor.

    Note that with a little planning, you will be in the situation such that if you re-source() your whole .R file, e.g., after correcting an error in the data or analysis, you will end up with a brand new, complete, readable report of the entire analysis with no effort.

    Note that the screen width affects the output by controlling the usual R text wrapping, e.g., with table(). Normally, you will want to keep the screen width around 60-70 characters to make it easier to read the report.

    Note that whenever you run report(x, new=TRUE, ...), if the report file name matches an existing file, the old file is deleted.

    Here are some examples that demonstrate what you can do:

         report("\nDemographics")
         report(table(age, gender))
         report(paste("Number of visits =", nrow(dat)))
         report("\nSuccess by treatment")
         report(with(dat, table(success, treatment, exclude=NULL)))
         report("\nYears of education:")
         report(summary(demog$educ))
         report(paste("\nDroppping", sum(noVisits|oneVisit), 
                      "subjects with no CERAD's or only 1 visit"))
         report(expression(str(my.data.frame)))
         

    The last example uses "expression()" because the "str()" function breaks the usual R rules and uses "cat()" rather than returning its result as an object. The "stem()" function is another example.

    There three helper functions in report.R.

    1. matForm(x, cols=12) converts a vector (string, numeric or factor) into a string matrix with a specific number of columns (even if length(x)%%cols!=0), so that long vectors don't ruin the appearance of the report.
    2. total(tab, margins=1:2) adds totals to the result of table()
    3. pct(tab, margins=1:2) adds percents to the result of table()

    Note that pct() and total() can both be used on the same table, in either order. Each respects the results of the other to avoid the incorrect and/or confusing output that could result from, e.g., including data and their total when computing percents.

    Optionally, you can use reportLatex() (code and description in reportLatex.R) to convert your .txt file into a .tex (Latex) file. This can incorporate plots as follows: when you are going through your analysis use the report text "See ... in myPlotFile.pdf", e.g.,

         plot(rnorm(20, type="b", main="Random normals", xlab="time", ylab="x")
         fname = "rnorm.pdf"
         dev.copy(pdf, fname); dev.off()
         # Important: Be sure to put a blank between "in" and the end quote
         #            since  sep="" will be in effect.
         report("\nSee 20 Gaussians in ", fname)
         

    With or without these special graphics commands, when you run reportLatex() a .tex file is created with the same base name as your .txt report file. Note that you can manually edit the .tex file at this point if desired.

    You then process this .tex file with pdflatex myReportFile.tex in Linux (or however else you know to process Latex files on any operating system) to produce the .pdf report file.

    If you used the special graphics indicator text "See ... in someFileName.pdf", then the plots will be included in the report, and the caption of the figure will be the text between "See" and "in". Also the caption will include figure numbers starting at "Figure 1".

    If you prefer to use a different graphics file type than "pdf" (as long as it is compatible with whatever version on pdflatex or latex that you are using) just run the optional form, e.g., reportLatex(extension=".pdf") substituting your graphics extension for "pdf".

Tips for Programming in R

1. End each function with return() or invisible() rather than using implicit returns. This conforms to standard programming practice in most other languages and make your program easier to read.

2. Start each function with checks of the arguments. It takes a little extra time but will usually repay you (or other users of the function) by pointing out the source of errors. Here is an example:

        myfun = function(dtf, name, p=0.5) {
          if (is.matrix(dtf)) dtf = data.frame(dtf)
          if (!is.data.frame(dtf)) stop("dtf must be a data.frame or matrix")
          if (!is.character(name) || length(name)!=1) stop("name must be a single character string")
          if (p<=0 || p>=1) stop("p must be in the interval (0,1)"
          ...
          return(rslt)
        }
        

3. Allow for stopping and restarting of functions with long loops (e.g., MCMC).
   A good trick is to setup your function (or even just a loop) as follows:

        myfun = function() {
          if (file.exists("myresults.dat")) {
             ...load and use old results...
          }
          ...
          for (i in 1:10000) {
            if (file.exists("stop")) {
              write.table(myresults, file="myresults.dat")
              stop("Early stop due to detection of stop file")
            }
            ...
          }
          ...
          return(...)
        }
        

   Then, you can create a file called "stop" at any time (e.g., in Linux using "echo stop>stop" at the Linux prompt) and the function will gracefully stop at the start of the next loop iteration. Without too much work, you can probably set up your function to automatically continue wherever you left off. Just remember to delete or rename the "stop" file before running the function again.

4. Avoid using "attach" as a way to save typing. The major problem is that modification of old elements or creation of new ones is not saved when you quit (and "save workspace") R. This leads to insidious errors. One alternative is to use with(), e.g., something like:

        with(mydtf, plot(x, y, col=gender))
        

   where the columns of "mydtf" are x, y, and gender.

5. Working with "non-visible functions": If you try, e.g., methods(logLik), you will find some methods (e.g., logLik.glm) that are marked with an asterisk and are "non-visible". Here is how to get a copy of those functions. Use getAnywhere(logLik.glm) to find that it is in the "namespace" called "stats". Then mylogLik.glm=stats:::logLik.glm will get you a copy of the function.

6. (Advanced) To make a nice user interface with dialog boxes, etc. consider the Tcl/Tk package. Here is a good introduction. Here is the R help. Here is a primer with an update. You might prefer a higher-level package called rpanel, described here and here, with this home page, and this package reference, and this cute little example which needs spacer.gif. Here are more R examples. Here are links about comparing tcl/tk to other systems. And here is a (non-R) Tcl/Tk Electronic Reference.

Problems and Solutions

1. Problem: After reading data from a file into a data.frame (e.g., with read.csv()) the first column name has a weird "ï.." prefix.     Solution: This is due to Windows adding a "byte order mark" to the beginning of the file. If you use, e.g., read.csv("myfile.csv", fileEncoding = "UTF-8-BOM"), the extra character will be removed, if necessary.
2. Problem: Loading dates, e.g., from Excel, and working with dates is poorly documented.     Solution: Load datetest.csv, then try the examples in Rdates.R.
3. Problem: Each click for the locator() annoyingly causes the bell to ring.
       Solution: options(locatorBell=FALSE)
4. Problem: Create a new data.frame column that is a complex code based on old columns.
       Solution: Create a function for one subject and apply() it to all subjects. This is much more efficient than a for loop. E.g.

        myfun = function(x) {
          # Argument x should contain one row, columns a,b,e,f.
          # The result is the mean of a and f unless b is missing or negative,
          # in which case the min of e and f is returned.
          if (is.na(x[2]) || x[2]<=0) {
            return(min(c(x[3],x[4])))
          } else {
            return((x[1]+x[4])/2)
          }
        }
        dtf$new = apply(dtf[,c("a","b","e","f")], 1, myfun)
        

   An alternative is as follows. The optional first line may prevent some wonky errors, and is good practice.

        dtf$new = NA  # in general, NA is safer than 0, protecting against bad logic
        Sel = is.na(dtf$b) | dtf$b<=0
        dtf[Sel, "new"] = pmin(dtf$e[Sel], dtf$f[Sel])
        Sel = !is.na(dtf$b) & dtf$b>0
        dtf[Sel, "new"] = (dtf$a[Sel]+dtf$f[Sel])/2
        

5. Problem: Analyze (all) subsets of a data.frame
      Solution: To analyze a single subset of a data.frame, you can use an index vector (logical or numeric) as the "row.selector" (first argument) of the form "incdata[row.selector, col.selector]". For example, the expression median(incdata[incdata$sex=="female", "income"]) calculates the median income of just the female subjects in data.frame "incdata".

   But expressions for each of several categories is awkward and inefficient. So the methods below present efficient alternatives. If no subsetting variable exists, consider using the R function cut() or Problem/Solution #2 to create it.

   Here is code you can paste into R to generate a sample data.frame to use as an example:

        n = 20
        incdata = data.frame(sex=c("male","female")[1+rbinom(n,1,0.5)],
                          race=c("black","white","hispanic","Asian")[1+rbinom(n,3,0.5)],
                          income=round(rnorm(n,50000,15000)),
                          networth=pmax(0,round(rnorm(n,50000,30000))))
        

   • The aggregate() command is designed to apply a built-in or user-defined function to each of one or several columns of a data.frame after automatically dividing the data.frame by the levels of one or more "by" variables.

          aggregate(incdata[,c("income","networth")], by=list(gender=incdata$sex,race=incdata$race), FUN=median)
          

     To get the results of a two-way aggregate into a table form, see aggregate.table() in the gdata package.

   • The tapply() function or a combination of the split() and lapply() functions can be used similarly to aggregate(). (In fact, aggregate() is a wrapper for tapply().) The split command can be applied to a vector, matrix, or data.frame, and will corresponding produce a list of vectors, matrices, or data.frames based on the splitting variable(s). Try, for example,

               with(incdata, split(income, race))
               split(incdata[,3:4], list(gender=incdata$sex, race=incdata$race))
               

     Two things to note that differ from using aggregate are that list() is not needed when a single split variable is used, and that the result includes empty elements for categories with no data. The second step is to use lapply() or sapply() on the split() list to carry out some function of the data, either a built-in function, a function that you have created, on an on-the-fly "anonymous" function. The difference between lapply() and sapply() is the latter will return a simple, non-list object, if possible. This example calculates the mean of the ratio of income to net worth for each sex/race subgroup using an anonymous function on the two-column elements of the list:

               tmp = split(incdata[,3:4], list(gender=incdata$sex, race=incdata$race))
               lapply(tmp, function(x) mean(x[,1]/x[,2]))
               sapply(tmp, function(x) mean(x[,1]/x[,2]))
               

     This example calculates range of income for each race:

               with(incdata, sapply(split(income, race), function(x) diff(range(x))))
               

     The simpler version uses tapply():

               with(incdata, tapply(income, race, function(x) diff(range(x))))
               

6. Problem: Need a function like read.table(), but the data are already in a character variable

       dat = c("category size type", "abc 23.4 g17a72", "aaa 3.2 h19h33", "bar 17 z12z12")
       dtf = read.table(textConnection(dat), header=TRUE)
     

7. Problem: R is too slow
   • First try some version of the "apply" command instead of "for" loops. (Note that sapply(1:n,FUN) can be used to run FUN n times if the first argument of FUN is a dummy (unused) argument.)
   • Consider using parallel computing, e.g., SNOW: (web page) (package) or Parallel R (paper) (package)
   • More experimental is Luke Tierney's R byte compiler. The help file is here. The zipped tar file is here. To load on Linux use something like gunzip foo.tar.gz to convert foo.tar.gz to foo.tar, then tar xvf foo.tar to extract all of the files. Use the Linux command R CMD INSTALL -l compiler compiler to install the package (-l compiler is needed only if you don't have admin permissions). In R, use library(compiler, lib.loc="compiler") to load the library for use (each time). (Change lib.loc= if the main directory is not named compiler in your current directory.)
   • Consider writing all or part of your code in C (see here).

Useful Functions

• Utility functions
  1. A function to Test Shiny Apps submitted as homework through the Canvas course management system. testShinyApps() modifies deals with the file renaming system of Canvas, including submission versions, and allows you to easily run and interact with Shiny apps that a class of students submits.
  2. Contrasts in R: SAS has easy to setup and use statements called CONTRAST and ESTIMATE for many PROC's. Something similar is needed for R. My function, contrast(), in contrast.R combines an lm() or glm() result with a contrast coefficient vector to produce a contrast estimate, standard error, confidence interval, and p-value. An example is at contrastExample.R. Note that named interactions can be entered using a form like L=c(facB=1, covar=-2, 'facB:covar=-2') when the model formula is Y~fac*covar.
  3. A string-to-value function that allows leading white space and trailing white space or non-numerics: strToVal() in strToVal.R
  4. hash.R is a function to mark output as comments by adding hash marks. This allows easy cut-and-paste into script files. E.g., hash(table(sample(LETTERS[1:5],10,rep=TRUE))) gives a result like:

       # A B C D E 
       # 3 2 2 2 1 
       

  5. Fix column names after read.csv() or read.table() to standard camel case using camelCase.R.
  6. R has cov2cor(), but the inverse function that computes the variance covariance matrix from a correlation matrix and a variance vector is not a part of standard R. Use cor2cov() in cor2cov.R.
  7. An R function to execute code pasted into the operating system clipboard is in cbr.R (ClipBoard Run). This can be used, e.g., to execute sample code found in a pdf file, whether or not it is commented or includes ">" and "+" prompt characters.
  8. An R function to save R output to the clipboard is in cbs.R (ClipBoard Save). Among other uses, this can be used to place R tables, including coefficient tables, into a Word or Excel document.
  9. An R function to save R tables to the clipboard is in dtfToTable.R. This is more full featured, including use of regular expressions to select rows and/or columns, combining of low/high CI values into a single CI expression of the form "[low, high]", and rounding with terminal zeros as needed.
  10. R code for SPSS versions of skewness and kurtosis and their standard errors: spssSkewKurtosis.R
  11. A function to convert tables to Latex is in table2Latex.R. This takes an R table and creates a file with the corresponding Latex code. The file can then be read into Latex using \input{filename}. The function has extensive support for custom control of bolding, cell justification, double lining, decimal places, and misc. adornments (e.g., color, italics, etc.).
  12. R is missing a function to extract the standard deviation and variance covariance matrix of an lme object. Use getRE() in getRE.R.
  13. Convert SPSS files to text files: Use getSPSS() in getSPSS.R.
  14. Convert various data files (including SPSS) to text files and/or a data.frame: Use convertData() in convertData.R.
  15. Package gdata has function read.xls() which can read Excel files. (Requires loading of PERL, e.g. ActivePerl.)
  16. Deprecated: Read Excel data into R using the ODBC method (with package RODBC): getExcel() in getExcel.R.
      Utility function getExcelStructure() reports sheet names and named ranges from an Excel workbook.
      Example: useExcel.R   test.xlsx
  17. Deprecated: Read Excel data into R using the COM method (with package rcom): dtfFromExcel() in dtfFromExcel.R.
  
  
• Pulling data from an image into an R data.frame If you see data in a plot image, e.g., in an article, and want to analyze it in R, use imageToData() from imageToData.R. To prepare the image for this function try these steps on a PC: Use PrintScreen to capture the image, load the image in the Paint program, crop the image, save as a 256-color bitmap. The use imageToData() function takes the name of the bitmap file and a pair of x and y coordinates marked on the image (to set the scale). See the header at the top of the file for more options. When you run the program, the image is displayed (but colors are not preserved) and you follow the instructions to click on the specified scale tic marks, and then on the points in the plot. The return data.frame contains the data. If the data have groups (e.g., different colors or symbols), you can use the function parameters groups= to specify either the number of group or their names. To try the example at the bottom of the file, use this example: imageToPlotExample.bmp (This function requires R package 'readbitmap'.)
• Genetics
  1. Convert a data frame column with entries like G/C to two columns with codes 1 and 2 with 0 for missing: SNPCodeSplit.R() in SNPCodeSplit.R
  2. Use MACH to extract the set of best haplotypes from HapMap (and optionally write input files for EHAP): bestHaplotypes() in BestHaplotypes.R
  3. Create input files for EHAP from a data.frame: writeEhap() in EhapRfunctions.R
  
  
• Simulation
  1. An interactive suite to learn about the Cox Proportional Hazards model through simulation: Home page for SimTool .
  2. A set of functions to simulate 1 or 2 way ANOVA data: ANOVAsim() in ANOVAsim.R.
  3. A set of functions to simulate ANOVA data for 1 between and 2 within subjects factors: ANOVAsim3() in ANOVAsim3.R.
  4. Simulate a mixture of four-parameter-logistic curves with random effects: sim4pl() in Sim4plRE.R.
  
  
• Exploratory Data Analysis and Residual Analysis
  1. A better interaction plotter: IAplot.R. This function works like interaction.plot(), but allows one to view the spread, either as individual points (e.g., to see outliers) or with 95% confidence intervals on the means. It is also very customizable.
  2. EDA for logistic regression: LREDA.R. This function plots the fraction of success vs. a continuous variable and optionally (using color) a categorical variable, with a smoothing curve. Counts (weights) are shown with different sized symbols. Note that the categorical variable could be a construct from two or more variables. Another version is BinContEDA() in BinContEDA.R. A version with confidence intervals is LREDACI.R.
  3. An interactive way to explore binned residual plotting vs. Hosmer Lemeshow goodness-of-fit for detecting non-linearity and interaction in logistic regression is implemented using the R package "shiny".
     Instructions:
     1. Start R in some directory
     2. Create a subdirectory called "LogReg"
     3. Place ui.R and server.R in the subdirectory
     4. Install package "shiny", if not yet installed
     5. Load shiny using library(shiny)
     6. Start the app with runApp("LogReg"). It will start in your browser (you may need to "permit" R to write to your browser).
     The app simulates data from a logistic model with two quantitative x values. You can add non-linearity by changing the "beta1^2" slider away from zero and/or changing the "beta12" interaction beta slide away from zero. The output includes:
     • a plot of the model fit for y~x1+x2 showing predicted probabilities vs. x1 at the 1/3 and 2/3 quantiles of x2 including 95% pointwise confidence bands (in black and red); on the plot the green curves show the probabilities for the true model
     • the binned residual plot (see Gelman and Hill's logistic regression chapter); radio buttons allow you to change the residual type and/or the number of bins
     • the Hosmer-Lemeshow statistic and p-value
     • the Hosmer-Lemeshow observed vs. fitted table
     • the model coefficients and standard errors
     As you change the beta's or other choices using the sliders or radio buttons on the left, the plots and other output on the right change automatically. To re-simulate data without changing the betas, click the "check box" called "Click to repeat" (each time you check or uncheck the box a new data set is simulated and displayed).
  4. A set of interactive functions to visually cluster latent class trajectory data: Home page for RVisClust.
  5. Residual plots for SAS PROC TRAJ: trajResiduals() and associated functions in trajResiduals.R
  6. Exploratory data analysis for lme(). Function lmeeda() in lmeeda.R plots slope vs. intercept for each subject in a groupedData object, optionally using different symbols for different levels of an additional variable.
  7. Extensive exploratory data analysis for "epochal" time series data: SurvEDA.R (test with SurvEDATest.R). This function produces color plots of several events per subject, with the option to line up on any one event, to sort by any variable or by an inter-event time, and to group by any variable.
  8. A 2-dimensional non-parametric confidence region plotter is xy95.R. This uses convex hulls to plot an approximate 2-d confidence region.
  9. A baseline hazard plotting function for R is baseHazPlot.R. Note that basehaz() produces the cumulative hazard, not the hazard.
  10. Weibull parameter estimation in a running window to look for secular change in parameter values: runWeibPlot() in RunningWeibulll.R
  11. Helper functions for package flexmix in the mixture of regressions context: FlexMixHelpers.R with examples in FlexSim.R. Includes trajectory plotting, re-ordering of latent classes, and computation and plotting of concomitant distributions.
  
  
• Other graphics
  1. Improved image plots (including a color legend and more control of axes): iplot() in iplot.R

  2. A simple, useful function for adding confidence intervals to plots is in addCIs.R.

  3. A full featured function to produce Gantt charts is in gantt.R.
  4. A function to make word clouds is in makeWordCloud.R
  5. Replacements for plot.factor() and barplot.default() to allow placement of legends for plot() with data in the form of factors. The original functions always place the legend at top right. The functions allow specification of lx= and/or ly= to place the legend (default is top right). Both arguments can be any value between 0 and 1. For lx, 0 means on the left, and 1 means on the right. For lx=0.5, the legend is centered; for lx<0.5 it is left justified; for lx>0.5 it is right justified. Vertical position is similarly specified using ly= with 1 being the top of the plot. FixPlotFactor.R

  6. A function to plot LDA or QDA borders: DAplot() in DAplot.R

  7. A function to produce means plots or 2- or 3-way interaction plots in R using the results from SAS GLIMMIX/MIXED/GLM is in SASIAPlot.R.

     This function handles arbitrarily complex models, and plots all (or some selected) combinations of any one, two, or three variables in the model, fixing all other variables at default or user selected values. The plotted means and 95% confidence intervals (actually, +/-2SE) are the same as you would get from writing (many) SAS ESTIMATE statements, but especially for models with higher order interaction and/or missing level combinations, this is extremely tedious in SAS. Note that many different useful plots can be produced from one SAS model run, by changing the R function's 'xVar', 'lineVar', and 'thirdVar' parameters, as well as the 'xLevels', 'lineLevels' and/or 'thirdLevels'

     This function produces useful plots just by specifying the two files containing the SAS output and the names(s) of the variables to plot. By setting additional function parameters, you have close control of the final plot details, including reordering, subsetting, and/or relabeling factor levels; relabeling the variable names; changing the title, and y-axis label; suppressing or keeping a subtitle which automatically shows all of the fixed values; changing font sizes, line colors, point characters, and line types, changing the location of the legends for the second (by color) and third (by line type) variables; keeping or suppressing the mean point markers; setting the width of the caps on the ends of the confidence bands, setting the y-axis value limits; and controlling the spread of the means of the levels of the second (line) and/or third variable around the nominal x-value for the first variable.

     For safety sake, the function writes fixed values that it chooses to the screen, as well as table(s) of levels vs. labels if the function thinks they may have been incorrectly set (i.e., labels do not appear similar to levels). The latter output can be suppressed with 'silent=TRUE'. Finally the whole plot can be suppressed with 'plot=FALSE' if you only want the return value of the means, standard errors, and upper and lower 95% confidence values.

     If the csv versions of the SAS ODS tables of coefficients and the variance-covariance matrix are named "coef.csv" and "covb.csv" you don't even need to supply those values to the function. And if no 'xVar', 'lineVar', and 'thirdVar' values are supplied, the first highest order plot available is produced by default.

     Categorical variables are fixed at their first values or user-specified values. Each quantitative covariate requires a fixed value in the 'fix' list or a set of values in the appropriate 'xLevels', 'lineLevels' and/or 'thirdLevels' parameters. As a convenience, 'xLabels', 'lineLabels' and/or 'thirdLabels' may be specified without the corresponding 'Levels' parameter if you just want to subset or reorder the levels, or even make simple (unambiguous) clearer labels.

     A full example, including the SAS code to generate the csv files is found in SASIAPlotTest.R. Although the original data needed to run the SAS code is not provided, the SAS output files needed to test the R function are supplied.

     Note that, because of a missing feature from PROC GLM, an extra step (creation of the ODS FIT table) is required. Therefore it may be easier to just fit the model with PROC MIXED without any random effects.

  8. A function to facilitate interactive use of R by providing clickable buttons and sliders is in buttonInput.R. See the examples at the end of the file.

     The basic idea is that you can flexibly specify a grid of buttons and sliders (possibly in groups) on some or all of a graphic window, then run a loop that receives information on what the user clicked so that you can execute any code you like, then loop back for more user input (repeatedly until the user uses the graph window 'Stop' button or your own stop button).
     The two basic modes of use are a menu in the margin of a plot for interactive plot control and a full screen menu with actions occurring in the console window (remember to use flush.console() in your code) or in one or more additional graphics windows.

  9. Weighted boxplots: wboxplot() has an argument weights= that is a vector of sample weights corresponding to the data to be plotted. For side-by-side boxplots, the weights should be a list of weight vectors. The argument scale= can be change (from a default of 1) to adjust the relative size of the outlier points. (Not all boxplots functionalities are currently included.) The function wboxplot() and related functions in WeightedBoxplot.R
     Or use Erich Neuwirth's wtd.boxplot() in package ENmisc from CRAN.

  10. A function to draw 95% empiric 2-dimensional confidence ellipsoids from sample data, e.g., a Bayesian posterior sample: xy95.R

  11. A function to draw 95% Gaussian 2-dimensional confidence ellipsoids from mean and covariance: ceGaussian2.R
  
  
• Power Calculation via Simulation
  • Power calculation for any model can be performed via simulation. Briefly, for any fixed set of conditions, you simulate data and compute the p-value many times and then report the power as the fraction of the time that p≤0.05.
  • The file powerSimClass.R contains convenience functions to implement running and reporting power calculations. An example is shown at the end of that file. These functions implement a new "S3" class called "powerSim" which is an object with power simulation results. The object is created with the function simPower() and takes as input a simulation function, an analysis function the returns a p-value, a range of sample sizes, a parameter vector, a coefficient index which points to which parameter element is being tested, a range of effect sizes for that parameter element, and a list of optional additional information.
  • Method functions are included so that print(), summary(), and plot() work directly on these objects. In addition, a merge() function merges the results from two powerSim objects into one new powerSim object, e.g., to fill in more sample sizes or effect sizes, or two improve precision by combining simulations from identical scenarios.
  • The functions may also be run as a "fake package". To do this, run attach("powerSimClass.RData") after placing powerSimClass.RData in an accessible location. This makes the functions available without having them appear in your workspace.

• WinBUGS/jags related
  1. My R package, rube, is a Really Useful WinBUGS Enhancer. It makes working with WinBUGS much easier. It is (currently) built on top of the R2WinBUGS package. R2WinBUGS takes a good first step to allow you to
     • Get your data and parameter initializations ready in R.
     • Run your model in WinBUGS without touching WinBUGS. (The doodle editor is really just a toy!)
     • Examine your results in R.

     The rube package is a wrapper for R2WinBUGS, but also much more. If you use it, you don't need to learn the details of R2WinBUGS.

     Note that when you install rube, it will install R2WinBUGS, if that is not already on your system. But you must install WinBUGS itself before using rube. Also note that you must upgrade R to at least version 2.10.0 before installing rube. With WinBUGS, the rube package only runs on a Windows operating system, because WinBUGS only runs on Windows. But if you use "jags", any OS will work. In that case, install the R package "rjags" and jags itself.

     To install rube see .

     Here are the main features of rube:

     1. A bugsCheck() function that does syntax checking of your model description "code" with explicit error messages. This finds most syntax errors and several other types of errors.
     2. While syntax checking only is the default for bugsCheck(), it can also help you find many additional problems related to the relationships between the model code, the data, and the initializations. When you set categorize=TRUE, a detailed summary of the above-mentioned relationships is produced.
     3. Easy to use "pseudo-functions" make your code much more flexible and easy to maintain.
        1. Anywhere in your code, you can use the "substitution" syntax VAR=value to set default values. This is most useful for hyperparameters. The two advantages are that it make the code easier to read (it becomes more obvious where you chose arbitrary hyperparameter values), and it makes hyperparameter sensitivity testing easy because you can add simple run-time options to override any or all defaults without modifying your code. The utility function showDefaults() can be used to list all of the defaults in a model. As an example,

           Y[i] ~ dnorm(0, Y.PREC=100)

           will be interpreted as

           Y[i] ~ dnorm(0, 100)

           if nothing special is done, but as

           Y[i] ~ dnorm(0, 1000)

           if subs=list(Y.PREC=1000) is used as a run-time function option.
        2. Coding of covariates is made simpler, easier to read, and much more flexible (while avoiding the maintenance problems of multiple model code files). This is accomplished through use of the linear combination pseudo-function which takes the form LC(prefix, FORMULA, suffix, index) anywhere in your code. Here prefix, suffix, and index are "hard-coded" in your model, but FORMULA is specified at run-time (with no default), and FORMULA is just like the right-hand-side of an R model formula (with some minor restrictions). The formula will include the covariates of interest for a particular WinBUGS model run, and the prefix and suffix are used to automatically create corresponding parameter names. An intercept is created unless -1 is included in the formula. Usually you will want to append the index to the covariates values. Here are a couple of examples:

           mu[j] <- dnorm(LC(b,DEMOGS,,j), Y.PREC=100)

           will be interpreted as

           mu[j] <- dnorm(beta0 + bOld*Old[j] + bMale*Male[j] + bOldMale*Old[j]*Male[j], 100)

           in various rube functions when you define varList=list(DEMOGS="Old*Male"). Then when you re-run your analysis with, say, varList=list(DEMOGS="Old*Male+White") the appropriate covariate model is automatically used. Your model can have any number of formulas, and the same name can be used in multiple places, or any number of different names can be used, depending on your model complexity. You can distinguish your parameters from the covariate data values using the prefix, suffix, or both.

           An example of a more complex covariate formula is "A+B+C+B:C+D*E*F+(G+H+I)^2" which is shorthand for the model formula

                            A + B + C + B:C + D + E + F + D:E + D:F + E:F + G + H + I + G:H + G:I + H:I
                            

           Currently the model formulas do not allow "-" (except -1 is allowed) or "/", and inside the (...)^n syntax, "A:B" and "A*B" are treated as "A+B".

           Warning: rube does not recode your variables, so the formula system will be inappropriate without warning for k-level categorical variables that have not been manually recoded into k-1 indicator variables (or other suitable coding). But the system does work for 0/1 or 1/-1 indicator variables and quantitative covariates.

        3. Parallel with LC() is FOR() which generates the corresponding prior distributions for the covariate parameters. The syntax is FOR(prefix, FORMULA, suffix, codeText). This syntax generates one model code line for each term/parameter in the FORMULA by replacing a question mark (?) in the "codeText" with the parameters one by one. Continuing with the above example, we might have a model code line

           FOR(b,DEMOGS,, ? ~ dnorm(0, COV.PREC=100))

           to automatically generate

                            b0 ~ dnorm(0, 100)
                            bOld ~ dnorm(0, 100)
                            bMale ~ dnorm(0, 100)
                            bOldMale ~ dnorm(0, 100)
                            bWhite ~ dnorm(0, 100)
                            

           when varList=list(DEMOGS="Old*Male+White") is defined. If intercept b0 has a different prior, code it separately and use varList=list(DEMOGS="Old*Male+White-1")
        4. A less used pseudo-function is BASE(FORMULA,index) which can be useful for encoding a baseline covariate status. E.g.

           isBase[i] <- BASE(DEMOGS,i)

           will be coded as

           isBase[i] <- (1-Old[i])*(1-Male[i])*(1-White[i])

           when varList=list(DEMOGS="Old*Male+White") is defined.
        5. The final improvement in the form of pseudo-functions is conditional coding. This allows two (or more) options two be in the code with the final choice made at run-time. One example of where this is useful is for a choice of distributions for a parameter. There are two forms of conditional coding: in-line and multi-line.

           In-line conditional coding takes the form code1 IFCASE(caseFormula) code2 ELSECASE code3 ENDCASE code4 which results in code1 code2 code4 if the caseFormula evaluates to TRUE and code1 code3 code4 if the caseFormula evaluates to FALSE. The code1, ELSECASE code3, and code4 elements are all optional. The "codeFormula" can be simple such as GAMMA or complex such as GAMMA || (ALPHA && !INTERCEPT) which will evaluate to FALSE and TRUE respectively when the argument subs=c("ALPHA","AGE2") is included in the various rube functions. This uses standard R evaluation including OR (||), AND (&&), NOT (!) and parentheses. Each element in the "subs" vector is defined TRUE and all missing elements are defined to be FALSE.

           The multi-line form is:

                             IFCASE(GAMMA || (ALPHA && !INTERCEPT))
                             a <- b + c
                             d <- e + f
                             ELSEIFCASE(ALPHA)
                             a <- b + c
                             d <- e - f
                             ELSECASE
                             a <- b
                             d <- f
                             ENDCASE
                             

           which evaluates to

                             a <- b + c
                             d <- e + f
                             

           when the argument subs=c("ALPHA","AGE2") is used, and

                             a <- b + c
                             d <- e - f
                             

           when the argument subs=c("ALPHA","INTERCEPT") is used.
     4. Rube provides a wrapper for the bugs() function of R2WinBUGS, called doBugs(), which includes the ability to generate data on the fly and to more flexible generate parameter starting values on the fly (see below), as well as extensive abilities to modify WinBUGS code on the fly (see above). Here is an example: Rasch.R.
     5. Rube is very flexible in the form of the WinBUGS model specification. It allows either a file with the specification or a string vector or a single vector separated by newline (\n) characters. The latter is convenient for keeping your model definition in the same file as your code, where you just define your model in this form:

                    myModel = "model {
                      for (i in 1:N) {
                      prec.y[i] <- pow(sig.y[i], -2)
                        Y[i] ~ dnorm(x[i], prec[i].y)
                      }
                    }"
                    

        To examine this very long string in R, use showModel(myModel). To see a particular version of the model based on evaluating any conditional code and/or pseudo-functions, use a form like

        showModel(myModel, subs=c("singleX","gamPrior"), varList=list(RHS1="A+B", RHS0="A+B-1"))

        Also note that showModel() with no arguments shows the most recently used form run through doBugs().
     6. Particularly for simulation studies, the ability to provide data as a data-generating function with optional arguments is a useful alternative to the usual method of providing data as a list.
     7. As is R2WinBUGS, parameter starting values can be set with either a data list or a function, but rube allows for passing arguments to this function. This allows alternate methods of starting value generation in different situations, possibly fully or partially data driven. In complex Bayesian problems, this can be invaluable.
     8. Rube provides a function called p3() which plots the trace, AR plot, and histogram for each monitored parameter, allowing backing up as well as quitting before seeing all parameters
     9. Rube provides a function called priPost() which is a prior/posterior density plotter to compare the posterior density of various distributions to their prior densities. This includes some "compound" prior distributions, e.g., a beta or gamma that has two gammas as its parameters.
  2. Functions for exploration of reasonable prior distributions are in hyperExplore.R (but will be moved to rube).
     • betaHyperExplore() allow interactive setting of the four hyperparameters of a beta distribution parameterized by two gamma distributions. A histogram of the prior density and an image plot of the means and sds of the betas that are being mixed to produce the final composite beta are shown. The final explored parameters are returned.
     • gammaHyperExplore() allow interactive setting of the four hyperparameters of a gamma distribution parameterized by two gamma distributions. A histogram of the prior density and an image plot of the means and sds of the betas that are being mixed to produce the final composite gamma are shown. The final explored parameters are returned.
     • betaHyperPlot() shows the histogram of a composite beta prior density from four gamma hyperparameters or the posterior histogram from the "bugs" object.
     • gammaHyperPlot() shows the histogram of a composite gamma prior density from four gamma hyperparameters or the posterior histogram from the "bugs" object.

All links active 5/15/2014. Please report missing links, errors, and suggestions to

To my Home Page