[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 http://www.nihrcrsu.org/ *t:*+44 (0)116 229 7268 *e:*ajs22 at le.ac.uk <mailto:ajs22 at le.ac.uk> *w:* www.le.ac.uk <http://www.le.ac.uk/>_ _cid:image001.gif at 01D0F6E3.BA906A20 Follow us on Twitter <https://twitter.com/uniofleicester> or visit our Facebook <http://www.facebook.com/uniofleicester> page__ cid:image002.png at 01D0F6E3.BA906A20 ** <http://www.avg.com/email-signature?utm_medium=email&utm_source=link&utm_campaign=sig-email&utm_content=emailclient&utm_term=oa-4885-b> Virus-free. www.avg.com <http://www.avg.com/email-signature?utm_medium=email&utm_source=link&utm_campaign=sig-email&utm_content=emailclient&utm_term=oa-4885-b> <#DAB4FAD8-2DD7-40BB-A1B8-4E2AA1F9FDF2>
_______________________________________________ R-sig-meta-analysis mailing list R-sig-meta-analysis at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-meta-analysis