Friday, June 25, 2010

debug() in R

While I'm on a (posting) roll, here is one more useful function, one that I'm a little embarrassed about discovering only now. It's debug() function in R. You tag you favorite function, foo(), with debug(foo). The next time foo() is executed, either by itself or is called from another function, debug executes foo() line by line, and provides an interactive prompt at each step so you can check values of the variables at the current evaluation, plot things, define things, or do whatever else you want. There are some built-in commands to be used at the prompt (if you can't think of something better):

n - advance to the nextstep
c - continue to the end of the current context: e.g. to the end of the loop if within a loop or to the end of the function
Q - quit debugging

In the end you can run undebug() to stop tracking the function.

For more, see here:

Taste of Chicago (off-topic)

Yeah, those of us who live here think it's lame, but I think I'm going to go anyway. They do have quite a few interesting restaurants (Ukrainian Shokolad Pastry & Cafe anyone?), so why not!

Here is a map (from Chicago Tribune):

fancyvrb package in LaTeX

In order to insert a code snippet (e.g. from R) into a LaTeX document, one could use a "verbatim" environment (e.g. with \begin{verbatim} and \end{verbatim}). This works, but isn't very customizable, e.g. cant' change font or fontsize easily.

Instead, download a package called fancyvrb. The manual is at  For example (and note capital V in Verbatim):


  function(X, tol = sqrt(.Machine$double.eps)) {
    ## Generalized Inverse of a Matrix
    dnx <- dimnames(X)
    if(is.null(dnx)) dnx <- vector("list", 2)
    s <- svd(X)
    nz <- s$d > tol * s$d[1]
      if(any(nz)) s$v[, nz] %*% (t(s$u[, nz])/s$d[nz]) else X,
      dimnames = dnx[2:1])


To install a LaTeX package (at least in my setup of WinEdt + MikTeX under Windows), do the following:
- Got to MikTeX Package Manager
- Wait for it to update the list of available packages
- Find the one you want to install, and double click on it.

Wednesday, June 23, 2010

Stata: stset and missing values

If a censoring indicator variable has missing values (but the actual time variable does not), one would expect Stata to omit that observation from subsequent survival analyses, or at least to give a warning that values are missing. However, it treats these observations as censored! For example: 

. tab ios, miss

        ios |      Freq.     Percent        Cum.
          0 |         53       46.09       46.09
          1 |         55       47.83       93.91
          . |          7        6.09      100.00
      Total |        115      100.00 

. count if osmo==. 

. stset osmo, fail(ios)

     failure event:  ios != 0 & ios < .
obs. time interval:  (0, osmo]
 exit on or before:  failure

      115  total obs.
        0  exclusions
      115  obs. remaining, representing
       55  failures in single record/single failure data
 2609.162  total analysis time at risk, at risk from t =         0
                             earliest observed entry t =         0
                                  last observed exit t =  93.92876 

. stci

         failure _d:  ios
   analysis time _t:  osmo

             |    no. of 
             |  subjects         50%     Std. Err.     [95% Conf. Interval]
       total |       115    29.42466      3.419791       22.126    35.9014

. tab  ios _d, miss

           |          _d
       ios |         0          1 |     Total
         0 |        53          0 |        53
         1 |         0         55 |        55
         . |         7          0 |         7
     Total |        60         55 |       115


Wednesday, June 16, 2010

Thursday, June 10, 2010

R: spaghetti plot using ggplot2 package

## color and group can be based on different variables


Friday, June 4, 2010

R: Spaghetti Plots for growth curves using lattice package

xyplot(vol~time, dataname, group = id, type="a")

xyplot(vol~time|trt, dataname, group=id, type="a")

SAS: Spaghetti Plots for growth curves using SAS v.9.2 PROC SGPLOT and PROC SGPANEL

Reposted from

proc sort data=dat1;
by trt id time;

proc sgplot data=dat1;
series x=time y=vol/group=trt break;

proc sgpanel data=dat1;
panelby trt;
series x=time y=vol/break group=trt;

UPDATE: how to embed pdf or ppt in the blog or webpage

If the PDF file that you want to embed is not on the web already and therefore does not have a URL that ends with *.pdf (e.g. it's a file on your computer), you can still embed it using third-party websites. I used Scribd, which is a document sharing site. 

Here are the steps:

1. Go to
2. Click the "Upload" button at the top of the screen
3. Enter your email (they promise not to spam)
4. Once your file is uploaded, click the "Share" link next to it, and copy the HTML link in the "Embed this document" section
5. Paste the link into your blogger post

Note, that this will create a Scrbd account for you and will email you password for future use. That's essentially it for how to create an account. 

Mathematica: Clearing values for typeset variables

I was surprised that a subscripted variable (e.g. Z_{k}) cannot be cleared using Clear[] function in Mathematica. Instead, you need to use Unset[]. See below.
Unset Subscripts                                                            

Thursday, June 3, 2010

how to embed pdf or ppt in the blog or webpage

Go to
and get the code

Stata 11.1

Stata 11.1 was released today (type update all to update, and then help whatsnew to see what's new). Of particular note given the past couple of posts on this blog is the following:

  • Mixed models. Linear mixed (multilevel) models have new covariance structures for the residuals: exponential, banded, and Toeplitz. See xtmixed.

lme4: Mixed-effects modeling with R

partial draft of the book by Douglas Bates

Older print version can be found here

Mixed-Effects Models in Sand S-PLUS
José C. Pinheiro and Douglas M. Bates

Subscribe via email

Enter your email address:

Delivered by FeedBurner


Blog Archive