## 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:

http://www.biostat.jhsph.edu/~rpeng/docs/R-debug-tools.pdf

### 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 http://ctan.org/tex-archive/macros/latex/contrib/fancyvrb/fancyvrb.pdf.  For example (and note capital V in Verbatim):

\begin{Verbatim}[fontsize=\small]

library(MASS)
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] structure( if(any(nz)) s$v[, nz] %*% (t(s$u[, nz])/s$d[nz]) else X,
dimnames = dnx[2:1])
}

\end{Verbatim}

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==.
0

. 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

### BBQ on Friday 7/30 after work (date changed again!)^3

location: Masha's new apartment

## Thursday, June 10, 2010

### R: spaghetti plot using ggplot2 package

## color and group can be based on different variables
qplot(time,vol,data=dataname,group=id,color=trt,geom="line")

qplot(time,vol,data=dataname,group=id,geom="line",facets=.~trt)

## 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

proc sort data=dat1;
by trt id time;
run;

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

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

### 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 http://www.scribd.com
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

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.

## Thursday, June 3, 2010

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-eﬀects modeling with R

partial draft of the book by Douglas Bates

http://lme4.r-forge.r-project.org/book/

Older print version can be found here

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

## Wednesday, June 2, 2010

### you might want to try -xtmixed- instead

As an alternative, the following link gives an excellent overview of some of the things that are possible with -xtmixed-: http://www.stata.com/meeting/fnasug08/gutierrez.pdf

### lmer for SAS PROC MIXED Users

http://cran.r-project.org/web/packages/SASmixed/vignettes/Usinglmer.pdf

LMER tutorial
http://www.stat.wisc.edu/~bates/PotsdamGLMM/LMMD.pdf
The following has R commands for plots
http://www.bioconductor.org/workshops/2008/PHSIntro/lme4Intro-handout-6.pdf