NYT: In Simulation Work, the Demand Is Real

The New York Times published this interesting article on how the ability to design and perform computer simulations is a highly marketable skill for careers across many disciplines.

In methodology development we use simulation nearly every day. We've developed our own specialized genetic data simulation software, genomeSIMLA, that's freely available here by request for PC, Mac, and Linux.

But if you have R on your computer (get it free here), here's how to do a really simple Monte Carlo simulation to determine the power of a one-sample t-test.

First, fire up R and type this command:

rnorm(100)

That command generates 100 random numbers drawn from a standard normal distribution, mean=0, sd=1. Now type this:

rnorm(100,mean=2,sd=7)

That also draws 100 random numbers from a normal distribution, but this time the mean is 2 and the standard deviation is 7. You can also get the same results by just typing this:

rnorm(100,2,7)

Now, let's do a one sample t-test:

t.test(rnorm(100,2,7))

That command performs a one-sample t-test on the 100 samples drawn from a normal distribution with mean=2 and sd=7. Remember, the null hypothesis of a one-sample t-test is usually "the mean is not significantly different from zero". So if the p-value is less than .05, we would typically reject this null, and say that the mean is significantly different from zero.

Now, we knew that the mean was different from zero, because we said draw from a distribution with mean=2. But if this was the case and we only drew 100 samples, how likely is it that we would detect a difference? That's the power of the test - given that the null is false, how likely is it that we reject the null hypothesis?

One way we can answer this question is with a simulation.

First, let's type the same command, but just get ONLY the p-value from the t-test:

t.test(rnorm(100,2,7))$p.value

Was is less than .05? Try typing it again (you can hit the up arrow key to bring up the last command in R, just like on the Linux command line). It will be different because we have a different set of 100 observations. Type it in over and over again. Sometimes it will be less than .05, other times it wont be. Let's do this 1000 times, and see how often it is less than .o5. Let's use the replicate command:

replicate(1000,t.test(rnorm(100,2,7))$p.value)

That simulates doing the t-test 1000 times, and gives you the p-value from each one.

Now, let's do a logical test to see which of those are less than .05:

replicate(1000,t.test(rnorm(100,2,7))$p.value)<0.05

If you typed that in you'll see lots of TRUE's and FALSE's. TRUE means that the t-test on that particular sample was less than .05. Now, internally, R represents TRUE as 1, and FALSE as 0. So if we take the average of all 1000 of these, that will tell us the proportion of times out of 1000 trials that the p-value of the one-sample t-test was less than .o5:

mean(replicate(1000,t.test(rnorm(100,2,7))$p.value)<0.05)

When I did this the power was right around 80%. If you do this again it will be slightly different because remember we are sampling randomly so the results will vary slightly!

Congratulations, you just did your first simulation / power study! Of course because we know what the null distribution of t-statistics looks like under the null, we can mathematically determine the power of a t-test without doing simulation studies:

power.t.test(n=100,delta=2,sd=7,sig.level=.05,type="one.sample")

But if we had developed our own method or algorithm we probably wouldn't have a mathematical formula to calculate power, which is why we rely on simulation studies. Be sure to check out my other posts on power calculation software, choosing the correct analyses, and code to run analyses in R and other software.

Side by side analyses in Stata, SPSS, SAS, and R

I've linked to UCLA's stat computing resources once before on a previous post about choosing the right analysis for the questions your asking and the data types you have. Here's another section of the same website that has code to run an identical analysis in all of these statistical packages, with examples to walk through (as they note here, just because they don't list code for a particular software for an analysis doesn't mean the software can't do it). They also have examples of power calculations for a smattering of statistics using previously mentioned G*power in addition to others.

UCLA Stat computing - data analysis examples with Stata, SPSS, SAS, R, and others

Get your genome sequenced by Illumina for $48k

This week Illumina launched their own personal genome sequencing service. For $48,000 they'll send you the sequence of your entire genome on a Mac computer that you can keep. According to their website, all the sequencing is done in a CLIA-certified clinical lab. One thing different about this than other consumer genetics services is that they require you to consult your doctor before signing up, and have them request sequencing for you, like writing a prescription. Then they send the sequence back to your doctor to discuss your results.

Now, even as a geneticist I'm not sure what to tell a layperson to do with all 3 billion of their bases sequenced, so what is a general practitioner to do when their patient seeks medical advice based on a service like this? Share your thoughts in the comments.

http://www.everygenome.com/

Challenges of translating genetic tests into clinical and public health practice

Here's a good paper published online this morning in Nature Reviews Genetics. As the title suggests, the paper covers many of the challenges associated with translating molecular genetic tests into clinical practice. There are tons of references to articles on the ethical, legal, and social issues surrounding genetic testing in health management. The authors also introduce and give reference to several methods from econometrics (value-of-information analysis, cost-benefit analysis, cost-effectiveness analysis, cost-utility analysis, etc) that are frequently used by policy and decision makers to prioritize technology and research funding. However, as the authors point out, although prioritizing funding of technology and research should ideally be based on rigorous evidence-based assessment and appraisal of the existing scientific evidence, this evidence is often not available or would otherwise be prohibitively expensive to collect. The authors conclude with a statement that prioritizing funding for translational genetic research remains just as much of a challenge as the genomic research itself.

Challenges of translating genetic tests into clinical and public health practice (NRG AOP June 9)

Make Pretty Regression Tables in Stata

The estout package for Stata is useful for quickly creating nicely formatted tables from a regression analysis for tables or papers. To install it, fire up Stata and type in this command:

ssc install estout, replace

Stata will automatically download and install the package. Run the regression as you normally would, then use the esttab command (part of the estout package) to create a table using those results. See their examples to see how the command works.

Click the thumbnail below to check out the difference. The top half is Stata's default output. The bottom is estout's formatted output.

Introductory statistics with R

I know that a lot of you are scrambling to spend your training grant money by next week. If you think you'll ever need to use R, I strongly recommend buying this book: Introductory Statistics with R, by Peter Dalgaard ($48, Amazon). I picked this up a while back and read through most of it in a day or two. It's not the best book to learn basic stats, but if you already have some background it's a great resource for learning how to do analyses in R, suitable for people who've never used R before. Incidentally, the author is an editor of the previously mentioned R journal.

You can actually view most of the content from an older edition of this book online for free at Google Books.

Use meaningful color schemes in your figures

Some of the best figure design ideas come from cartographers. If you've ever read a Tufte book you've seen lots of examples. Let's talk about using color effectively. Penn State geography professor Cindy Brewer's ColorBrewer tool for selecting color schemes for figures has been conveniently packaged into an R library called RColorBrewer. You'll have to read up on the RColorBrewer documentation to see how it works, but here I just want to point out how to use color schemes in a meaningful way based on the type of data you're presenting.

In general, three useful color schemes are diverging, sequential, and qualitative. The example palettes shown here can be reproduced with the following R code:

install.packages("RColorBrewer")
library(RColorBrewer)
display.brewer.all(type="seq")
display.brewer.all(type="div")
display.brewer.all(type="qual")

A diverging color scheme is useful for de-emphasizing the mean value, or for drawing attention to departures from a critical midpoint in either direction, such as in a normal distribution. Here's an example color palette and an example of how you might use a diverging scheme:





A sequential color scheme is useful for de-emphasizing the zero, or the lower bound to the data, while highlighting the importance of increasing values. Here's a palette and an example barchart.





Finally, a qualitative color scheme provides high contrast between adjacent values, and is useful for categorical or nominal data. Here's a palette and for an example check out the picture from the previously mentioned genetic diversity in Africa paper:




See Cindy Brewer's explanation of this topic for more details and examples of each scheme.

A good GWAS review for the uninitiated

GWAS reviews are a dime a dozen these days, but I found this one in particular that's a good up-to-date review suitable for people relatively new to the field. While this one has a focus on psychiatric genetics, it has a short summary on topics like positional methods, candidate gene studies, common variation, rare variation, CNVs, GWAS design, and issues concerning power and sample size. There's a table defining some basic terminology, a timeline of methods from linkage to the 1000 genomes project, and a table summarizing potentially problematic GWAS design issues.

Genomewide Association Studies: History, Rationale, and Prospects for Psychiatric Disorders (Pubmed)

For other links to send new folks joining your group, check out posts tagged with Tutorials or Recommended Reading.