[R-meta] Reproducing results using regtest in metafor
James Pustejovsky
jepu@to @end|ng |rom gm@||@com
Fri Aug 9 20:04:49 CEST 2019
Option 1 amounts to fitting the model
Y_i = b0 + b1 m.c_i + b2 se_i + e_i
with the assumption that Var(e_i) = sigma^2 v_i for unknown sigma^2. It
can be fit using the nlme package using the following syntax (building off
of the previous code):
library(nlme)
gls_fit <- gls(yi ~ m.c + sqrt(vi), data = ma.dataset, weights = varFixed(~
vi))
summary(gls_fit)
I think the correct df would be 27 in this case.
With option 2, the residuals from the first stage have lost 2 degrees of
freedom (down to 28), and then they lose two more with the second stage fit
(because the intercept is re-estimated when it should be constrained to
zero).
James
On Fri, Aug 9, 2019 at 8:24 AM Sutton, Alex (Prof.) <ajs22 using leicester.ac.uk>
wrote:
> Dear Michael
>
> Many thanks for taking the time to write. (And let me thank you for all
> your hard work on the R meta-analysis web page - this is a fantastic
> resource I use all the time)
>
> I think your suggestion is very plausible and I was thinking along related
> lines. Then I convinced myself that the df used in first regression was
> reflected in the se of the residuals (in a very "hand wavy" sort of way!) .
>
> As you say, the issue then focuses on which is "better" - my concern with
> option 1 is that I don't know how to reproduce it outside of the package -
> but as you say it may well be theoretically superior.
>
> I am also hoping Wolfgang can supply some definitive insight!
>
> Thanks again
>
> Alex
>
> Alex Sutton
> Professor of Medical Statistics
>
> Department of Health Sciences
> College of Life Sciences
> University of Leicester
> George Davies Centre
> 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 using le.ac.uk
> w: www.le.ac.uk
>
> Follow us on Twitter or visit our Facebook page
>
>
>
>
> -----Original Message-----
> From: Michael Dewey [mailto:lists using dewey.myzen.co.uk]
> Sent: 09 August 2019 13:19
> To: Sutton, Alex (Prof.) <ajs22 using leicester.ac.uk>;
> r-sig-meta-analysis using r-project.org
> Cc: DOLEMAN, Brett (UNIVERSITY HOSPITALS OF DERBY AND BURTON NHS
> FOUNDATION TRUST) <brett.doleman using nhs.net>; Freeman, Suzanne C. (Dr.) <
> suzanne.freeman using leicester.ac.uk>
> Subject: Re: [R-meta] Reproducing results using regtest in metafor
>
> Dear Alex
>
> I suspect the answer is that regtest.rma refits the model and so "knows"
> you have a moderator and hence loses you a degree of freedom whereas
> regtest.default does not "know" that and so proceeds as though you had just
> the single predictor (sei). That probably explains the change in t too
> although we may need to wait for Wolfgang to give us the definitive answer.
>
> Of course that leaves open the question which is better. I think your
> option 1 is what I would prefer but that is just based on instainct rather
> than any higher mathematics.
>
> Michael
>
> On 08/08/2019 17:42, Sutton, Alex (Prof.) wrote:
> > Hi All list members
> >
> > I wish to adjust a meta-analytic dataset by the control arm event rate
> > using meta-regression (ignoring any regression to the mean issues) and
> > then run Egger's funnel asymmetry test after the adjustment.
> >
> > (I paste complete data, R code and output below my query.)
> >
> > The regtest command in metafor can do this and an example is given in
> > the manual (option 1 in the code).
> >
> > But a colleague tried to do the same analysis in another package -
> > STATA
> > - but did not get the same result.
> >
> > I have tried to reconcile the 2, and have succeeded getting the STATA
> > result in R using the code in option 2. But I don't understand why
> > option 1 and option 2 code below give different results - I thought
> > they should be equivalent? (The degrees of freedom are different as
> > well as the test p-values)
> >
> > The key lines extracted from the below are:
> >
> > model2 <- rma(yi, vi, data = ma.dataset, mods = cbind(m.c),
> > control=list(maxiter=200, stepadj=0.5))
> >
> > print(summary(model2))
> >
> > # Option 1: Put model 2 directly into regtest
> >
> > option1 <- regtest(model2, model = "lm", predictor="sei")
> >
> > print(option1)
> >
> > # Option 2: Extract residuals and SE from model 2 and then put these
> > into regtest manually
> >
> > a <- rstandard.rma.uni(model2)
> >
> > option2 <- regtest(x=a$resid, sei=a$se, model = "lm", predictor="sei")
> >
> > print(option2)
> >
> > Hoping someone can offer me some insight.
> >
> > Many thanks in advance.
> >
> > Alex
> >
> > DATA BELOW:
> >
> > n.c
> >
> >
> >
> > n.t
> >
> >
> >
> > m.c
> >
> >
> >
> > m.t
> >
> >
> >
> > sd.c
> >
> >
> >
> > sd.t
> >
> > 15
> >
> >
> >
> > 15
> >
> >
> >
> > 43.34736
> >
> >
> >
> > 22.6778
> >
> >
> >
> > 20.14634
> >
> >
> >
> > 13.2187
> >
> > 15
> >
> >
> >
> > 15
> >
> >
> >
> > 36.66456
> >
> >
> >
> > 26.95877
> >
> >
> >
> > 14.24828
> >
> >
> >
> > 14.89395
> >
> > 15
> >
> >
> >
> > 15
> >
> >
> >
> > 51.74782
> >
> >
> >
> > 28.44658
> >
> >
> >
> > 23.83604
> >
> >
> >
> > 26.70313
> >
> > 15
> >
> >
> >
> > 15
> >
> >
> >
> > 49.80753
> >
> >
> >
> > 19.04398
> >
> >
> >
> > 24.16307
> >
> >
> >
> > 11.40324
> >
> > 15
> >
> >
> >
> > 15
> >
> >
> >
> > 60.67109
> >
> >
> >
> > 44.62121
> >
> >
> >
> > 29.01649
> >
> >
> >
> > 23.60531
> >
> > 15
> >
> >
> >
> > 15
> >
> >
> >
> > 31.34018
> >
> >
> >
> > 14.65917
> >
> >
> >
> > 18.27904
> >
> >
> >
> > 13.30416
> >
> > 25
> >
> >
> >
> > 25
> >
> >
> >
> > 44.48585
> >
> >
> >
> > 22.18729
> >
> >
> >
> > 24.98121
> >
> >
> >
> > 11.94804
> >
> > 25
> >
> >
> >
> > 25
> >
> >
> >
> > 42.19917
> >
> >
> >
> > 27.97304
> >
> >
> >
> > 27.69334
> >
> >
> >
> > 19.4292
> >
> > 25
> >
> >
> >
> > 25
> >
> >
> >
> > 44.99957
> >
> >
> >
> > 26.77364
> >
> >
> >
> > 33.80657
> >
> >
> >
> > 18.46828
> >
> > 25
> >
> >
> >
> > 25
> >
> >
> >
> > 62.6953
> >
> >
> >
> > 36.15569
> >
> >
> >
> > 32.1401
> >
> >
> >
> > 28.75357
> >
> > 25
> >
> >
> >
> > 25
> >
> >
> >
> > 35.38889
> >
> >
> >
> > 19.77613
> >
> >
> >
> > 16.96604
> >
> >
> >
> > 15.69643
> >
> > 25
> >
> >
> >
> > 25
> >
> >
> >
> > 36.46263
> >
> >
> >
> > 27.69418
> >
> >
> >
> > 22.18271
> >
> >
> >
> > 22.04955
> >
> > 25
> >
> >
> >
> > 25
> >
> >
> >
> > 42.63761
> >
> >
> >
> > 26.83291
> >
> >
> >
> > 22.60301
> >
> >
> >
> > 23.47935
> >
> > 50
> >
> >
> >
> > 50
> >
> >
> >
> > 57.92533
> >
> >
> >
> > 30.96353
> >
> >
> >
> > 28.90083
> >
> >
> >
> > 17.40233
> >
> > 50
> >
> >
> >
> > 50
> >
> >
> >
> > 26.15427
> >
> >
> >
> > 21.06599
> >
> >
> >
> > 17.39962
> >
> >
> >
> > 14.29371
> >
> > 50
> >
> >
> >
> > 50
> >
> >
> >
> > 39.98736
> >
> >
> >
> > 29.13268
> >
> >
> >
> > 19.99763
> >
> >
> >
> > 18.44498
> >
> > 50
> >
> >
> >
> > 50
> >
> >
> >
> > 46.11657
> >
> >
> >
> > 27.34703
> >
> >
> >
> > 24.23706
> >
> >
> >
> > 15.55458
> >
> > 50
> >
> >
> >
> > 50
> >
> >
> >
> > 48.0687
> >
> >
> >
> > 34.85292
> >
> >
> >
> > 25.19593
> >
> >
> >
> > 21.33037
> >
> > 50
> >
> >
> >
> > 50
> >
> >
> >
> > 53.43902
> >
> >
> >
> > 31.4177
> >
> >
> >
> > 31.33719
> >
> >
> >
> > 19.16536
> >
> > 50
> >
> >
> >
> > 50
> >
> >
> >
> > 56.2192
> >
> >
> >
> > 39.338
> >
> >
> >
> > 28.10021
> >
> >
> >
> > 24.66106
> >
> > 50
> >
> >
> >
> > 50
> >
> >
> >
> > 29.25929
> >
> >
> >
> > 15.60332
> >
> >
> >
> > 17.00783
> >
> >
> >
> > 12.57914
> >
> > 50
> >
> >
> >
> > 50
> >
> >
> >
> > 30.72385
> >
> >
> >
> > 20.22988
> >
> >
> >
> > 18.93253
> >
> >
> >
> > 15.11565
> >
> > 15
> >
> >
> >
> > 15
> >
> >
> >
> > 55.86973
> >
> >
> >
> > 27.38925
> >
> >
> >
> > 33.32485
> >
> >
> >
> > 22.29828
> >
> > 15
> >
> >
> >
> > 15
> >
> >
> >
> > 47.10929
> >
> >
> >
> > 28.39797
> >
> >
> >
> > 21.93013
> >
> >
> >
> > 25.12405
> >
> > 15
> >
> >
> >
> > 15
> >
> >
> >
> > 31.87787
> >
> >
> >
> > 15.31057
> >
> >
> >
> > 17.05947
> >
> >
> >
> > 10.99065
> >
> > 15
> >
> >
> >
> > 15
> >
> >
> >
> > 36.91457
> >
> >
> >
> > 26.14574
> >
> >
> >
> > 18.46971
> >
> >
> >
> > 18.57387
> >
> > 15
> >
> >
> >
> > 15
> >
> >
> >
> > 49.84609
> >
> >
> >
> > 34.04816
> >
> >
> >
> > 24.64221
> >
> >
> >
> > 21.22332
> >
> > 25
> >
> >
> >
> > 25
> >
> >
> >
> > 57.06665
> >
> >
> >
> > 32.91816
> >
> >
> >
> > 27.6288
> >
> >
> >
> > 19.84755
> >
> > 25
> >
> >
> >
> > 25
> >
> >
> >
> > 50.57364
> >
> >
> >
> > 36.87816
> >
> >
> >
> > 33.87053
> >
> >
> >
> > 21.38415
> >
> > 25
> >
> >
> >
> > 25
> >
> >
> >
> > 28.80046
> >
> >
> >
> > 21.97428
> >
> >
> >
> > 16.45233
> >
> >
> >
> > 15.75454
> >
> > CODE BELOW:
> >
> > library(metafor)
> >
> > rm(list=ls())
> >
> > ##simdat <- read.csv("s4_data.csv")
> >
> > ma.dataset <- escalc(n1i = n.t, n2i = n.c, m1i = m.t, m2i = m.c,
> >
> > sd1i = sd.t, sd2i = sd.c, data = simdat, measure
> > = "MD",
> >
> > append = TRUE)
> >
> > model2 <- rma(yi, vi, data = ma.dataset, mods = cbind(m.c),
> > control=list(maxiter=200, stepadj=0.5))
> >
> > print(summary(model2))
> >
> > # Option 1: Put model 2 directly into regtest
> >
> > option1 <- regtest(model2, model = "lm", predictor="sei")
> >
> > print(option1)
> >
> > # Option 2: Extract residuals and SE from model 2 and then put these
> > directly into regtest
> >
> > a <- rstandard.rma.uni(model2)
> >
> > option2 <- regtest(x=a$resid, sei=a$se, model = "lm", predictor="sei")
> >
> > print(option2)
> >
> > OUTPUT BELOW:
> >
> > Mixed-Effects Model (k = 30; tau^2 estimator: REML)
> >
> > logLik deviance AIC BIC AICc
> >
> > -84.1061 168.2122 174.2122 178.2088 175.2122
> >
> > tau^2 (estimated amount of residual heterogeneity): 0 (SE =
> > 6.3793)
> >
> > tau (square root of estimated tau^2 value): 0
> >
> > I^2 (residual heterogeneity / unaccounted variability): 0.00%
> >
> > H^2 (unaccounted variability / sampling variability): 1.00
> >
> > R^2 (amount of heterogeneity accounted for): 100.00%
> >
> > Test for Residual Heterogeneity:
> >
> > QE(df = 28) = 17.1497, p-val = 0.9456
> >
> > Test of Moderators (coefficient(s) 2):
> >
> > QM(df = 1) = 24.9726, p-val < .0001
> >
> > Model Results:
> >
> > estimate se zval pval ci.lb ci.ub
> >
> > intrcpt 3.5641 3.8211 0.9327 0.3510 -3.9251 11.0534
> >
> > m.c -0.4633 0.0927 -4.9973 <.0001 -0.6450 -0.2816 ***
> >
> > ---
> >
> > Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> >
> >>
> >
> >>
> >
> >> # Option 1: Put model 2 directly into regtest
> >
> >> option1 <- regtest(model2, model = "lm", predictor="sei")
> >
> >> print(option1)
> >
> > Regression Test for Funnel Plot Asymmetry
> >
> > model: weighted regression with multiplicative dispersion
> >
> > predictor: standard error
> >
> > test for funnel plot asymmetry: t = -0.2787, df = 27, p = 0.7826
> >
> >> # Option 2: Extract residuals and SE from model 2 and then put these
> >> directly into regtest
> >
> >> a <- rstandard.rma.uni(model2)
> >
> >> residuals <- a$resid
> >
> >> SEresiduals <- a$se
> >
> >>
> >
> >> option2 <- regtest(x=residuals, sei=SEresiduals, model = "lm",
> >> predictor="sei")
> >
> >> print(option2)
> >
> > Regression Test for Funnel Plot Asymmetry
> >
> > model: weighted regression with multiplicative dispersion
> >
> > predictor: standard error
> >
> > test for funnel plot asymmetry: t = -0.2170, df = 28, p = 0.8297
> >
> > *Alex Sutton
> > Professor of Medical Statistics
> >
> > ***
> >
> > Department of Health Sciences
> >
> > College of Life Sciences
> >
> > University of Leicester
> >
> > George Davies Centre
> >
> > University Road
> >
> > LEICESTER LE1 7RH UK
> >
> > *Please use Lancaster Road, Leicester, LE1 7HA for SatNav*
> >
> > Member of the Complex Reviews Support Unit
> >
> > https://eur03.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.n
> > ihrcrsu.org%2F&data=02%7C01%7Cajs22%40leicester.ac.uk%7Ce117781ca2
> > a44ca568b608d71cc3bd96%7Caebecd6a31d44b0195ce8274afe853d9%7C0%7C1%7C63
> > 7009499341114715&sdata=0iAykupucHNp6qBZHrBdJiw9BLww5gTGlTsWDIFg2Zo
> > %3D&reserved=0
> >
> > *t:*+44 (0)116 229 7268
> > *e:*ajs22 using le.ac.uk <mailto:ajs22 using le.ac.uk>
> > *w:*
> > https://eur03.safelinks.protection.outlook.com/?url=www.le.ac.uk&d
> > ata=02%7C01%7Cajs22%40leicester.ac.uk%7Ce117781ca2a44ca568b608d71cc3bd
> > 96%7Caebecd6a31d44b0195ce8274afe853d9%7C0%7C1%7C637009499341114715&
> > ;sdata=i4CJas8E39vy%2BcRWETeeFLVYCndvYsojfoKxwVwLT3c%3D&reserved=0
> > <https://eur03.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.
> > le.ac.uk%2F&data=02%7C01%7Cajs22%40leicester.ac.uk%7Ce117781ca2a44
> > ca568b608d71cc3bd96%7Caebecd6a31d44b0195ce8274afe853d9%7C0%7C1%7C63700
> > 9499341114715&sdata=D02HtRTW%2BRjXayhrM6tkUlvvOY21cPt2FT4Yo5qfzFE%
> > 3D&reserved=0>_
> > _cid:image001.gif using 01D0F6E3.BA906A20
> > Follow us on Twitter
> > <https://eur03.safelinks.protection.outlook.com/?url=https%3A%2F%2Ftwi
> > tter.com%2Funiofleicester&data=02%7C01%7Cajs22%40leicester.ac.uk%7
> > Ce117781ca2a44ca568b608d71cc3bd96%7Caebecd6a31d44b0195ce8274afe853d9%7
> > C0%7C0%7C637009499341114715&sdata=mxn3i9TylT7fLGKv44ULs3xQjYpxULu2
> > 1X6guJqxiuY%3D&reserved=0> or visit our Facebook
> > <https://eur03.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.
> > facebook.com%2Funiofleicester&data=02%7C01%7Cajs22%40leicester.ac.
> > uk%7Ce117781ca2a44ca568b608d71cc3bd96%7Caebecd6a31d44b0195ce8274afe853
> > d9%7C0%7C0%7C637009499341114715&sdata=JakhmIjpCe2VXycfx9ZYxdMT1z4r
> > Kx6wv1fyrfzdJpg%3D&reserved=0> page__
> >
> >
> >
> > cid:image002.png using 01D0F6E3.BA906A20
> >
> > **
> >
> >
> > <
> https://eur03.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.avg.com%2Femail-signature%3Futm_medium%3Demail%26utm_source%3Dlink%26utm_campaign%3Dsig-email%26utm_content%3Demailclient%26utm_term%3Doa-4885-b&data=02%7C01%7Cajs22%40leicester.ac.uk%7Ce117781ca2a44ca568b608d71cc3bd96%7Caebecd6a31d44b0195ce8274afe853d9%7C0%7C1%7C637009499341114715&sdata=G94bbFUIO8t%2FguFLoO1HFHodjMF3mVGLFiM7ejuFIko%3D&reserved=0>
>
> > Virus-free.
> > https://eur03.safelinks.protection.outlook.com/?url=www.avg.com&da
> > ta=02%7C01%7Cajs22%40leicester.ac.uk%7Ce117781ca2a44ca568b608d71cc3bd9
> > 6%7Caebecd6a31d44b0195ce8274afe853d9%7C0%7C1%7C637009499341124710&
> > sdata=aYlIcQAyezpanV387TaeLwlIqki0jpJ0%2FDspYgWNvfk%3D&reserved=0
> > <https://eur03.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.
> > avg.com%2Femail-signature%3Futm_medium%3Demail%26utm_source%3Dlink%26u
> > tm_campaign%3Dsig-email%26utm_content%3Demailclient%26utm_term%3Doa-48
> > 85-b&data=02%7C01%7Cajs22%40leicester.ac.uk%7Ce117781ca2a44ca568b6
> > 08d71cc3bd96%7Caebecd6a31d44b0195ce8274afe853d9%7C0%7C1%7C637009499341
> > 124710&sdata=cnQ7G%2Fd9EBWdWuKkagE9enGenkla%2FCpUF4X6T3IaXxU%3D&am
> > p;reserved=0>
> >
> >
> > <#DAB4FAD8-2DD7-40BB-A1B8-4E2AA1F9FDF2>
> >
> > _______________________________________________
> > R-sig-meta-analysis mailing list
> > R-sig-meta-analysis using r-project.org
> > https://eur03.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat
> > .ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-meta-analysis&data=02%7C01%7
> > Cajs22%40leicester.ac.uk%7Ce117781ca2a44ca568b608d71cc3bd96%7Caebecd6a
> > 31d44b0195ce8274afe853d9%7C0%7C1%7C637009499341124710&sdata=uZQWN8
> > g5lR%2Flox9rs2YExt4Czj%2B8lRl9Ch0j563nAyE%3D&reserved=0
> >
>
> --
> Michael
>
> https://eur03.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.dewey.myzen.co.uk%2Fhome.html&data=02%7C01%7Cajs22%40leicester.ac.uk%7Ce117781ca2a44ca568b608d71cc3bd96%7Caebecd6a31d44b0195ce8274afe853d9%7C0%7C1%7C637009499341124710&sdata=eIw%2FeMVcYsc8gQBrJkDF77GUTZMZxehWneBmAPG%2FhpQ%3D&reserved=0
>
> _______________________________________________
> R-sig-meta-analysis mailing list
> R-sig-meta-analysis using r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-meta-analysis
>
[[alternative HTML version deleted]]
More information about the R-sig-meta-analysis
mailing list