Meta analysis in R

or 

the beneficial effect of teaching on research

I have been fascinated by meta analysis for a long time. It is so obviously the right way to approach the true effect of an intervention. Recently, an old binder presented itself in a pile of shite bunch papers I meant to read but found myself throwing in the bin tidying away. It contained the draft of a database of physiological data from the first years of my PhD. The idea was to compare all the baseline data from our kidney research group in Uppsala to look at the effect of the models as such and the interventions that were used repeatedly. With the insight of the intervening years it seems a lot less interesting now, but I still have the feeling that some areas of experimental research could benefit from meta analysis.

Which brings us to the the story I am about to tell. Four or five years ago, when I moved back to Uppsala, I got offered a lecture on physiological changes in the elderly. It was to be part of a final year course for the master in biomedicine programme. Just a single hour to show how all the physiology from the rest of the programme changed with age. To compare and contrast ageing as such with the accrued ailments of living for a long time, and distinguish these from the chronic and age related disease. It was not a huge success, but given the title I was not too disillusioned. The second year I was given two hours. Still a bit on the short side, but one hundred percent better than one.

It is not the most popular lecture, but I have had it for five years now and one of the things I teach is that some parts of ageing is caused by metabolism itself. The burning of oxygen singes the organism and with time it will break much like the paneling in an old sauna. As proof of this I used the idea of caloric restriction, which can prolong life in many strains of yeast, mice, and rats. Then, in 2012 an article was published on the effect of caloric restriction in the Rhesus monkey, a primate, and reasonably the closest relative to humans in which an experiment could be expected to be finished any time soon. It showed no effect. I happily included this in my lecture as a counter-point. Until in 2014, when updating the lecture for a new semester, I found that another experiment with caloric restriction in Rhesus monkeys had published their data and found a clear difference.

This made it hard to continue the lecture as I had done, I could just show both studies and say that we don't know. But the total number of animals included was quite large, and the effect measure very straight-forward. Death. So, I performed a meta analysis of mortality of these two studies, and a third smaller study published in 2003. This is the story of that analysis.

Quickly I installed the R package rmeta by Thomas Lumley and set to work. It is quite easy really, we start with setting up a table of results from the included studies. The table should include the total number of subjects in each group, and the number of deaths per group.

Hultström, M. Acta Physiol (Oxf). 2015 Feb 14. doi: 10.1111/apha.12468.

The we push this trough the rmeta function meta.MH(). To get a forest plot, we just run the plot() command, which has a default for handling the result of meta.MH() in the form of a forest plot. If you have a larger meta analysis there is also the funnelplot() that can be used to assess publication bias. Anyway, the result is quick and easily understood, which is really one of the major strengths of the forest plot.

Hultström, M. Acta Physiol (Oxf). 2015 Feb 14. doi: 10.1111/apha.12468.

There was no significant effect of caloric restriction on all cause mortality in Rhesus monkeys. Or, rather there was a small, clearly non-significant, effect. One of the reviewers asked what would be needed to show if this effect was true. That is, could I please perform a power analysis. So, I installed the pwr package and ran a 2p2n.test() using the most generous effect estimate, i.e. a hypothetical study that ran to completion where the whole control population had died giving an effect of 0.08. This resulted in a required population of 2806 subjects to reach 85% power. This is the power-level which is normally used as the basis for power calculation in clinical studies. However, the age-related mortality was a different story that you can find in the actual article.

The next thing that surprised me was how difficult it was to get this simple little analysis published. It appears that experimental journals don't publish meta analyses, and clinical journals that publish meta analyses, don't publish experimental results. Finally, I found a benevolent editor at Acta Physiologica who permitted it to be published as an editorial. So that is where it resides today, and finally I can give a fairly clear answer in my lecture on the effect of reducing metabolism by caloric restriction on ageing and on mortality. Only problem is, I now have to explain meta analysis and forest plots before I can show the actual data.

And, no I am not going to starve myself so that I can avoid some diseases we can treat in favour for a frailty for which the only known treatment is eating more.

Winter-meeting in Oslo

The Norwegian Hypertension Society held its bi-annual scientific meeting in Oslo last week. It was a long time in planning. We set the place and date about two years ago. Sent the first announcement in the spring of 2012 and the call for abstracts in September. Then, at the very deadline, the abstracts and registrations start to trickle in. As with all meetings, we pushed the submission-deadline back a week just to allow people without basic planning skills with pressed schedules to join the meeting and present their data. Then, in as little time as possible, all the abstracts have to be formatted into a program, and session-chairs has to be found and matched so that they have an interest but aren't speaking themselves.

Anyway, we got some 17 free communications, varying from experimental physiology to international research politics, but with a heavy focus on epidemiology and clinical research. There were a couple of talks on the recently very hot topic of renal sympathetic denervation for treatment-resistant hypertension (more on that in another post), some interesting sub-group analyses from the LIFE and SCAST studies, and follow-ups on the now 40-year-old Oslo-Ischemia-Study. More of the program at the society home-page.

In addition, we had two invited lectures on statistics in clinical research. The first on how to develop and validate prognostic models by Ingar Holme, and the second on over-adjustment bias in multiple regression models held by Knut Liestøl. It was a useful repetition of the uses and pitfalls of these two very similar kinds of models that require very different study-designs and give very different information in the end. A common problem is that one tries to get etiological information from prognostic research, i.e. treating a risk-factor as a cause for the chosen end-point even though the observational design makes that impossible. The converse is equally common, i.e. trying to infer prognostic information from etiological studies, such as clinical trials, where the highly selected population makes general conclusions very suspect.

All told, it was a very successful meeting, well worth the time both for planning it, and attending.

Speech impediment - Locale in R under OSX

My R-installation developed a language problem ("L" is a statistical ploglamming language). The startup message, all errors and system messages started appearing in German. I do speak German but it was still quite annoying. This ended with a work-around that I published previously in a post about R-startup scripts. I have finally figured out what causes the error. My Language & Text settings for language were as follows:

1. British English 
2. Swedish
3. Norwegian
4. Danish
5. German 
6. French 

The problem lies in using the British English locale. I use it because I hate when all programs always correct my excellent spelling from Proper English™ to American. R does not have a locale for British English, nor for Swedish, Norwegian or Danish. So, to R, the locale priority list actually looks like this:

1.  
2. 
   
3. 
   
4. 
   
5. German 
6. French

So, understandably it goes for what it knows. The solution is to extend the list and put a standard English locale high up:

1. British English
2. English 
3. Swedish
4. Norwegian
5. Danish
6. German 
7. French 

Analysis of genome-scale data

I thought I should do a post on strategies for handling genome-scale data from large experiments. Like many happy novices we rushed in and performed a couple of huge microarray experiments. Now we are trying to compare these against eachother and against strain-differences in in the form of SNPs. To say that this is not trivial is an understatement. In the end we have come up with a couple of strategies that might be useful to others as well. They might seem simple, but they will help you wade through your data, without getting stuck there for too long.

Don't look at untested data for anything but quality control. In most analysis programs it is very easy to just plot a gene by its intensity, either over all samples or by crude difference between group-averages. This is fun, and can give you some idea of what you are looking for. However, there are anything from 20k-30k genes in your dataset and you will not have time to look at all of them and make a reasoned decision. That means that if you pick genes this way you will indeed pick some hits, but not in a systematic way, and that leaves you open to errors. Either of picking a gene where the variation is too large, or that may only be the n+100th gene in your set and not the top one. The point is: You just don't know. The suggestion is to test your genes first, and then use the list of significantly differentially expressed genes to pick your top hits.

Consider doing a single comparison across several groups. If you have run a large experiment, with several groups (where several is anything from four to whatever). It is a chore to keep track of which genes were significant in which comparison, and even more to then pick a small number of genes to validate that are somehow representative across groups. So, try to pool your data by phenotype. If you have one group where something happens and lots of different controls. Just pool the controls and do a single comparison. That will give you just one list and make your selection process much easier.

When comparing experiments, use rankings. So, now you have a list, or several if pooling didn't make sense. Now you have to find your top candidates for validation. It is then helpful to use rankings. The easiest (and often quite accurate) way is to rank by average difference between groups (that would be your fold-change or log-ratio). If you have just one list, then you are done. If you have two, or several lists, then you can order your lists by rank-sum or rank-product. The rank-sum is a kind method that gives you some hits even if a gene is poorly-ranked (high number) in one of your groups. The rank product will pronounce those that are highly-ranked in all lists. Then you have the very particular rank-difference and rank-ratio. These may be very interesting if you have two experiments with wildly different phenotypes, but that you can not compare directly. If a gene is ranked in the top 10 in one experiment and 7815th in the other, it may be a very interesting gene and will then be ranked highly by the rank-difference or ratio.