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 builtin 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/Rdebugtools.pdf
Useful biostatistics and statistical programming tips, comments, discussions, opinions, etc.
Friday, June 25, 2010
Taste of Chicago (offtopic)
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):
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/texarchive/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.
Instead, download a package called fancyvrb. The manual is at http://ctan.org/texarchive/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
. 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
Friday, June 18, 2010
Wednesday, June 16, 2010
BBQ on Friday 7/30 after work (date changed again!)^3
location: Masha's new apartment
comment and RSVP below, please
comment and RSVP below, please
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)
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~timetrt, dataname, group=id, type="a")
xyplot(vol~timetrt, dataname, group=id, type="a")
SAS: Spaghetti Plots for growth curves using SAS v.9.2 PROC SGPLOT and PROC SGPANEL
Reposted from http://statcode.blogspot.com/2009/12/sasspaghettiplotsforgrowthcurves.html
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;
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 thirdparty 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
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.
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
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.
Thursday, June 3, 2010
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: Mixedeﬀects modeling with R
partial draft of the book by Douglas Bates
http://lme4.rforge.rproject.org/book/
Older print version can be found here
MixedEffects Models in Sand SPLUS
José C. Pinheiro and Douglas M. Bates
http://www.springerlink.com/content/9781441903174/
http://lme4.rforge.rproject.org/book/
Older print version can be found here
MixedEffects Models in Sand SPLUS
José C. Pinheiro and Douglas M. Bates
http://www.springerlink.com/content/9781441903174/
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
Subscribe to:
Posts (Atom)
Blog Archive

▼
2010
(66)

▼
June
(16)
 debug() in R
 Taste of Chicago (offtopic)
 fancyvrb package in LaTeX
 Stata: stset and missing values
 BBQ on Friday 7/30 after work (date changed again!...
 R: spaghetti plot using ggplot2 package
 R: Spaghetti Plots for growth curves using lattice...
 SAS: Spaghetti Plots for growth curves using SAS v...
 UPDATE: how to embed pdf or ppt in the blog or web...
 Mathematica: Clearing values for typeset variables...
 how to embed pdf or ppt in the blog or webpage
 Stata 11.1
 lme4: Mixedeﬀects modeling with R
 you might want to try xtmixed instead
 lmer for SAS PROC MIXED Users

▼
June
(16)