[R-meta] Reproducing results using regtest in metafor

Fri Aug 9 14:18:53 CEST 2019

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 using le.ac.uk <mailto:ajs22 using le.ac.uk>
> *w:* www.le.ac.uk <http://www.le.ac.uk/>_
> _cid:image001.gif using 01D0F6E3.BA906A20
> Follow us on Twitter <https://twitter.com/uniofleicester> or visit our 
> Facebook <http://www.facebook.com/uniofleicester> page__
> 
> 	
> 
> cid:image002.png using 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 using r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-meta-analysis
> 

-- 
Michael
http://www.dewey.myzen.co.uk/home.html