ひさかたぶりの投稿
こういうデータをメタアナリシスで解析する。
testが処理した群、contが対照群
testで病気になった人がtest_posで
全体がtest_total
id test_pos test_total cont_pos cont_total 1 test1 5 177 6 144 2 test2 5 132 4 155 3 test3 2 23 33 47 4 test4 6 76 22 86 5 test5 10 146 2 134 6 test6 1 32 2 56 7 test7 2 133 4 155 8 test8 4 59 12 66 9 test9 1 48 1 44 10 test10 4 104 2 87
通常は、
meta1 <- metabin(test_pos, test_total, cont_pos, cont_total, data = dat, sm = "OR", studlab = id, comb.fixed = T, comb.random = T) meta1 OR 95%-CI %W(fixed) %W(random) test1 0.6686 [0.1998; 2.2372] 9.26 12.00 test2 1.4862 [0.3908; 5.6521] 5.10 11.32 test3 0.0404 [0.0083; 0.1960] 28.51 10.08 test4 0.2494 [0.0951; 0.6540] 27.38 13.29 test5 4.8529 [1.0436; 22.5681] 2.80 10.29 test6 0.8710 [0.0759; 10.0001] 2.03 6.57 test7 0.5763 [0.1039; 3.1975] 5.24 9.44 test8 0.3273 [0.0993; 1.0782] 15.21 12.08 test9 0.9149 [0.0555; 15.0827] 1.47 5.52 test10 1.7000 [0.3039; 9.5111] 3.02 9.40 Number of studies combined: k=10 OR 95%-CI z p-value Fixed effect model 0.5156 [0.3471; 0.7658] -3.2815 0.001 Random effects model 0.6065 [0.2703; 1.3607] -1.2129 0.2252 Quantifying heterogeneity: tau^2 = 1.0368; H = 1.69 [1.21; 2.37]; I^2 = 65% [31.2%; 82.2%] Test of heterogeneity: Q d.f. p-value 25.69 9 0.0023
で図にすると
となる。
metaforパッケージだと
metafor1 <- escalc(measure="OR", ai = test_pos, bi = test_total - test_pos, ci = cont_pos, di = cont_total- cont_pos, data = dat, append = T) metafor1res <- rma(metafor1$yi, metafor1$vi, data=dat, method ="DL") forest(metafor1res) Random-Effects Model (k = 10; tau^2 estimator: DL) tau^2 (estimated amount of total heterogeneity): 1.0348 (SE = 0.7940) tau (square root of estimated tau^2 value): 1.0173 I^2 (total heterogeneity / total variability): 64.92% H^2 (total variability / sampling variability): 2.85 Test for Heterogeneity: Q(df = 9) = 25.6562, p-val = 0.0023 Model Results: estimate se zval pval ci.lb ci.ub -0.5002 0.4121 -1.2138 0.2248 -1.3078 0.3074
となり、図にすると
となる。
で、気になるのは、このheterogenityとかが何を指すのか。
だけど、普通の混合モデルをlme4でやってみるために、まずはREMLでmeta analysisをやってみる。
こうなる。
metafor2res <- rma(metafor1$yi, metafor1$vi, data=dat, method ="REML") > metafor2res Random-Effects Model (k = 10; tau^2 estimator: REML) tau^2 (estimated amount of total heterogeneity): 1.1297 (SE = 0.8331) tau (square root of estimated tau^2 value): 1.0629 I^2 (total heterogeneity / total variability): 66.89% H^2 (total variability / sampling variability): 3.02 Test for Heterogeneity: Q(df = 9) = 25.6562, p-val = 0.0023 Model Results: estimate se zval pval ci.lb ci.ub -0.4971 0.4240 -1.1724 0.2410 -1.3281 0.3339
となる。
で、lme4でやってみると、
id syori pos total test1 test 5 177 test1 cont 6 144 test2 test 5 132 test2 cont 4 155 test3 test 2 23 test3 cont 33 47 test4 test 6 76 test4 cont 22 86 test5 test 10 146 test5 cont 2 134 test6 test 1 32 test6 cont 2 56 test7 test 2 133 test7 cont 4 155 test8 test 4 59 test8 cont 12 66 test9 test 1 48 test9 cont 1 44 test10 test 4 104 test10 cont 2 87
とデータをいじってから
lme1 <- glmer(cbind(pos, total) ~ syori + (1|id), data = dat2, family = binomial) Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod'] Family: binomial ( logit ) Formula: cbind(pos, total) ~ syori + (1 | id) Data: dat2 AIC BIC logLik deviance df.resid 122.5 125.5 -58.2 116.5 17 Scaled residuals: Min 1Q Median 3Q Max -1.9940 -0.6087 -0.2307 0.6815 2.2630 Random effects: Groups Name Variance Std.Dev. id (Intercept) 0.9123 0.9552 Number of obs: 20, groups: id, 10 Fixed effects: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.6772 0.3327 -8.048 8.42e-16 *** syoritest -0.5222 0.2027 -2.576 0.01 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Correlation of Fixed Effects: (Intr) syoritest -0.217
試験ごと、処理ごとにランダム効果としてみても
> lme2 <- glmer(cbind(pos, total) ~ 1+ (1|id) +(1|syori), data = dat2, family = binomial) > lme2 Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod'] Family: binomial ( logit ) Formula: cbind(pos, total) ~ 1 + (1 | id) + (1 | syori) Data: dat2 AIC BIC logLik deviance df.resid 125.7663 128.7535 -59.8831 119.7663 17 Random effects: Groups Name Std.Dev. id (Intercept) 0.9762 syori (Intercept) 0.2849 Number of obs: 20, groups: id, 10; syori, 2 Fixed Effects: (Intercept) -2.93
似たような数値はない。
調べてみると、
http://www.metafor-project.org/doku.php/tips:rma_vs_lm_and_lme
どうやらメタアナリシスはオッズ比なりリスク比を元にランダム効果モデルを当てはめているようですね。
???うーん？BUGSで式書いてやってみないと整理できないな。
