Wolfgang Viechtbauer, Ph.D., Statistician | Department of Psychiatry and
Neuropsychology | Maastricht University | P.O. Box 616 (VIJV1) | 6200 MD
Maastricht, The Netherlands | +31 (43) 388-4170 | http://www.wvbauer.com
-----Original Message-----
From: R-sig-meta-analysis [mailto:r-sig-meta-analysis-bounces at r-project.org] On Behalf Of Sutton, Alex (Prof.)
Sent: Wednesday, August 23, 2017 18:27
To: r-sig-meta-analysis at r-project.org
Cc: Sutton, Alex (Prof.)
Subject: [R-meta] Reproducing an R meta-analysis in STATA
Dear Group (I appreciate I maybe stretching outside the remit of this group a little with this request)
I have done an analysis in R for a publication, but the (reasonably stats savvy) collaborating clinician wants to know how to do the same analysis in STATA (which he uses and he thinks most of the readers of the paper would use). I am stuck on one specific point and wondered if anyone has the knowledge to help me.
The analysis adjusts a meta-analysis dataset for a moderator variable (using a mixed effects model) and then plots a funnel plot of the (raw) residuals vs. the standard error of the residuals.
The "metafor" package in R makes this very easy using code of the form below:
res <- rma(yi, vi, data = ma.dataset, mods = cbind(moderator))
funnel(res, level=0)
Now, in STATA, there is a "metareg" macro to fit the same mixed effects model, but no subsequent command to automatically plot the funnel plot of the residuals vs se(residuals). So I need to extract / compute these quantities manually (then there are funnel type commands - no problem). The raw residuals themselves are no problem to calculate (predicted values come as default from the "predict" command and then subtract them from the observed data) but I have spent a (embarrassingly) long time failing to compute the se(residuals)!
I have to confess not being able to find the formulae for se(residuals) in the literature anywhere (can anyone advise?) but in examining the source code to "metafor" believe a distilled version of the relevant code is:
ImH <- diag(x$k) - H
ve <- ImH %*% tcrossprod(x$M,ImH)
sei <- sqrt(diag(ve))
Which I understand to be
ve = (I-H)x M x transpose(I-H),
where M is a matrix with diagonal equal to the variance estimates of each study and 0 otherwise, H is the hat matrix and I is the identity matrix.
Then se(residuals) ?= the square root of the diagonal elements of the "ve" matrix.
Now I don't think I can get STATA to predict se(resid) for me via the predict command so I went about this via matrix commands to create matrices and compute the above (as is done in the R code). I can create everything but the Hat matrix. Predict can give the diagonal elements of this but not the whole matrix(!).
So I guess (on a good day) I could (possibly) work out how to construct the hat matrix myself, but the idea was to publish the simplest code possible to encourage non-statisticians to use it and this is just going in the opposite direction.
Am I being dense / missing something, are there alternatives to using the hat matrix to calculate se(resid)? Or is this a case where it simply isn't possible in STATA without more programming?
Many thanks in advance for any advice you can give
Best wishes
Alex
Alex Sutton
Professor of Medical Statistics
Department of Health Sciences,
College of Medicine, Biological Sciences and Psychology,
University of Leicester, Centre for Medicine, University Road, Leicester, LE1 7RH, UK
Please use Lancaster Road, Leicester, LE1 7HA for SatNav
Member of the Complex Reviews Support Unit
http://www.nihrcrsu.org/
t: +44 (0)116 229 7268
e:? ajs22 at le.ac.uk
w: www.le.ac.uk
Follow us on Twitter or visit our Facebook page