[R-meta] Computing Effect Size for Difference in Differences with Different Populations

Reza Norouzian rnorouz|@n @end|ng |rom gm@||@com
Fri Mar 17 05:03:48 CET 2023


Dear Mika,

You can certainly use both effect size measures and see to what extent the
choice can affect the results.

But if the traditions in your area of research allow, it may also be worth
running your model using the usual standardized mean differences (SMDs) as
the effect size of choice for each group at each time point and then
conduct appropriate contrasts to get the difference in differences (DID)
meta-analytically.

I believe in the **dev. version** of metafor, Wolfgang has recently added a
new function called emmprep() that in conjunction with library emmeans may
be helpful to achieve this relatively easily.

For example, in the data below (consisting of 3 groups and 3 time points),
you can compute the DID across 3 groups for long-term gains (i.e.,
differences in gains between time1 and time3) as follows:

# install dev version of metafor
library(emmeans)

dat$effect <- 1:nrow(dat)

# Model:
a1 <- rma.mv(smd ~group*time, V=vi, random = list(~time|study, ~1|effect),
data=dat,
            dfs = "contain")

# New function in metafor:
sav1 <- emmprep(a1)

# Marginal means:
emms <- emmeans(sav1, ~group*time)

# Desired contrasts among marginal means:
contrast(emms, list("Long-term: Gain(gr2) - Gain(gr1)" =
c(1,-1,0,0,0,0,-1,1,0),  # length must match #of rows in emms
                                "Long-term: Gain(gr2) - Gain(gr3)" =
c(0,-1,1,0,0,0,0,1,-1),
                                 "Long-term: Gain(gr3) - Gain(gr1)" =
c(1,0,-1,0,0,0,-1,0,1)
             ))

The above steps assume that you are interested in comparing the groups at
the final stage of the studies (time3), but see differences in their
initial status (time1) that demand taking those differences into account.
Similar steps can be taken to target differences at time2 taking into
account the initial differences.

The above procedure could be used with other models as well. For instance,
pretending that dat is not multilevel in structure, we could fit the
following model and repeat the above steps.

a2 <- rma(smd ~group*time, vi=vi, data=dat, test = "knha")

Kind regards,
Reza

# Data
dat <- structure(list(study = c(1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3,
                                    3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5,
5, 5, 5, 5, 5, 5, 5, 5,
                                    5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 5, 6, 6, 6, 6, 6, 6,
                                    7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 9,
9, 10, 10, 10, 10, 10,10),
                          smd = c(1.681, 4.122, 2.6, 3.457, 1.546, 3.072,
0.22, 3.545,
                                4.454, 0.624, 0.047, 0.306, 1.48, 1.2,
1.483, 1.054, 0.638, 0.998,
                                -0.31, -0.079, 0.955, 0.721, 1.531, 0.974,
-0.104, 0.235, -0.11,
                                -0.005, -0.117, 0.384, -0.194, 0.299,
0.452, 0.66, 0.268, 0.132,
                                0.078, 0.254, -0.339, -0.057, 0.482, 0.278,
0.104, 0.253, 0.107,
                                0.652, 0.386, 0.134, -0.329, -0.366,
-0.011, -0.359, 0, -1.279,
                                -0.178, -0.066, -0.149, -0.547, -0.173,
-0.279, 0.795, 0.172,
                                0.422, 0.461, 1.591, 1.001, -0.352, 0.473,
0.453, 0.3, 0.798,
                                0.868, 0.897, 0.862),
                           vi = c(0.208, 0.481, 0.284, 0.384, 0.2,
                                  0.335, 0.08, 0.206, 0.278, 0.114, 0.106,
0.101, 0.138, 0.125,
                                  0.127, 0.124, 0.111, 0.112, 0.072, 0.073,
0.08, 0.077, 0.092,
                                  0.081, 0.048, 0.049, 0.048, 0.048, 0.047,
0.047, 0.047, 0.047,
                                  0.049, 0.051, 0.049, 0.048, 0.048, 0.049,
0.049, 0.048, 0.048,
                                  0.047, 0.047, 0.047, 0.047, 0.049, 0.047,
0.047, 0.112, 0.112,
                                  0.11, 0.112, 0.11, 0.132, 0.101, 0.099,
0.101, 0.103, 0.101,
                                  0.1, 0.196, 0.168, 0.186, 0.172, 0.239,
0.189, 0.136, 0.138,
                                  0.149, 0.135, 0.157, 0.146, 0.16, 0.146),
                          group = structure(c(1L, 2L, 1L, 2L, 1L, 2L, 1L,
1L, 1L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L,
                                              3L, 1L, 2L, 1L, 2L, 1L, 2L,
1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 1L,
                                              1L, 1L, 1L, 1L, 1L, 1L, 1L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L,
                                              1L, 1L, 1L, 1L, 1L, 1L, 2L,
1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,
                                              2L, 1L, 1L, 1L, 2L, 1L, 2L,
1L, 2L), levels = c("1", "2", "3"
                                              ), class = "factor"),
                         time = structure(c(1L, 1L, 2L, 2L, 3L, 3L, 1L, 2L,
3L, 1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 1L, 1L, 2L,
                                            2L, 3L, 3L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 2L, 2L, 2L, 2L, 3L,
                                            3L, 3L, 3L, 2L, 2L, 2L, 2L, 3L,
3L, 3L, 3L, 1L, 1L, 2L, 2L, 3L,
                                            3L, 1L, 1L, 2L, 2L, 3L, 3L, 1L,
1L, 2L, 2L, 3L, 3L, 1L, 2L, 1L,
                                            1L, 2L, 2L, 3L, 3L), levels =
c("1", "2", "3"), class = "factor")),
                         class = "data.frame", row.names = c(NA, -74L))


On Thu, Mar 16, 2023 at 5:47 PM Mika Manninen via R-sig-meta-analysis <
r-sig-meta-analysis using r-project.org> wrote:

> Hey Wolfgang,
>
> Thank you very much for the reply.
>
> My mistake with the sd1i argument, I wasn't supposed to be using the post
> sds. Thank you for pointing that out.
>
> In my data the pre-treatment SDs are not similar (neither are the
> pre-treatment means) between the groups. That is why I am a bit unsure
> regarding the most appropriate ES computation method. I have mostly
> meta-analysed RCT data in which the two groups are almost identical at
> pre-treatment. I could not really find examples in which the difference in
> treatment response is being compared between two different populations
> receiving the same treatment.
>
> In any case, the following is an okay representation of the actual data and
> depending on which ES computation approach I use, the result looks quite
> different. So, I was wondering if you could help me in deciding which of
> the two is more appropriate or perhaps there is a third option. The data
> (or a subsection of it) can also be meta-analysed with raw effect sizes
> which does lead to different conclusions as well compared to the SMCRH/SMCC
> approach).
>
> Thank you very much in advance,
> Mika
>
> ### dataset
> set.seed(123)
> n_G1 <- rpois(50, lambda = 50)
> n_G2 <- n_G1
>
> postm_G1 <- rnorm(50, mean = 14, sd = 2.5)
> prem_G1 <- rnorm(50, mean = 10, sd = 2)
> postsd_G1 <- rnorm(50, mean = 2.4, sd = 0.4)
> presd_G1 <- rnorm(50, mean = 1.9, sd = 0.3)
>
> postm_G2 <- rnorm(50, mean = postm_G1 - 7.5, sd = 1.8)
> prem_G2 <- rnorm(50, mean = prem_G1 - 5, sd = 1.2)
> postsd_G2 <- rnorm(50, mean = 1.2, sd = 0.2)
> presd_G2 <- rnorm(50, mean = 1, sd = 0.2)
>
> G <- data.frame(prem_G1,presd_G1, postm_G1, postsd_G1, n_G1,
> prem_G2,presd_G2, postm_G2,postsd_G2, n_G2)
> G
>
>
> # Option 1 (could be SMCC as well I suppose)
> G1 <- escalc(measure="SMCRH", m1i=postm_G1, m2i=prem_G1, sd1i=presd_G1,
> ni=n_G1, sd2i = postsd_G1, ri=c(rep(0.7,50)), data=G)
> G2 <- escalc(measure="SMCRH", m1i=postm_G2, m2i=prem_G2, sd1i=presd_G2,
> ni=n_G2, sd2i = postsd_G2, ri=c(rep(0.7,50)), data=G)
> dat <- data.frame(yi = G1$yi - G2$yi, vi = G1$vi + G2$vi)
> dat
> # Crude mean of Effect sizes
> mean(dat$yi)
>
> # Option 2
> pldpre_sd = sqrt((presd_G1^2 + presd_G2^2) / 2)
> ES = ((postm_G1 - prem_G1) - (postm_G2 - prem_G2)) / pldpre_sd
> ES
> # Crude mean of Effect sizes from this formula
> mean(ES)
>
>
> to 16. maalisk. 2023 klo 17.36 Viechtbauer, Wolfgang (NP) (
> wolfgang.viechtbauer using maastrichtuniversity.nl) kirjoitti:
>
> > Hi Mika,
> >
> > Depends on what you mean by 'best'. But note that escalc(measure="SMCRH,
> > ...) computes (m1i-m2i)/sd1i, so you are using the post-treatment SD to
> > standardize, which is a bit unusual and different from your second
> approach
> > where you use the average pre-treatment SD to standardize. This aside,
> the
> > two approaches should lead to rather similar estimates, especially if the
> > pre-treatment SDs (assuming those are used in both approaches) are
> similar
> > across the two groups.
> >
> > Best,
> > Wolfgang
> >
> > >-----Original Message-----
> > >From: R-sig-meta-analysis [mailto:
> > r-sig-meta-analysis-bounces using r-project.org] On
> > >Behalf Of Mika Manninen via R-sig-meta-analysis
> > >Sent: Monday, 13 March, 2023 17:28
> > >To: R meta
> > >Cc: Mika Manninen
> > >Subject: [R-meta] Computing Effect Size for Difference in Differences
> with
> > >Different Populations
> > >
> > >Dear community,
> > >
> > >I am currently working on a meta-analysis that aims to examine the
> > >difference in training effects between two populations. Both
> > >populations underwent the same training, but at pre-test, the groups
> > >have significantly different means and standard deviations (about
> > >1-2sd difference in means).
> > >
> > >I am interested in computing the effect size for the difference in
> > >differences between the two groups. Specifically, I would like to know
> > >what is the best way to calculate the effect size given the
> > >significant difference in means and standard deviations at pre-test.
> > >
> > >Would the below be roughly accurate (Option 1):
> > >
> > >Option 1.
> > >
> > >G1 <- escalc(measure="SMCRH", m1i=postm_G1, m2i=prem_G1,
> > >sd1i=postsd_G1,ni=n_G1, sd2i = presd_G1, ri=c(rep(0.7,10)), data=G)
> > >G2 <- escalc(measure="SMCRH", m1i=postm_G2, m2i=prem_G2,
> > >sd1i=postsd_G2, ni=n_G2, sd2i = presd_G2, ri=c(rep(0.7,10)), data=G)
> > >dat <- data.frame(yi = G1$yi - G2$yi, vi = G1$vi + G2$vi)
> > >
> > >Option 2.
> > >
> > >ES = (G1 post_mean - G2 pre_mean) - (G2 post_mean - G2 pre_mean) /
> > pldpre_sd
> > >pldpre_sd = sqrt((presdG1^2 + presdG2^2) / 2)
> > >
> > >### dataset
> > >
> > >set.seed(123)
> > >
> > >postm_G1 <- rnorm(100, mean = 14, sd = 2.5)
> > >prem_G1 <- rnorm(100, mean = 10, sd = 2)
> > >postsd_G1 <- rnorm(100, mean = 1.4, sd = 0.2)
> > >presd_G1 <- rnorm(100, mean = 1, sd = 0.2)
> > >n_G1 <- rpois(100, lambda = 50)
> > >
> > >postm_G2 <- rnorm(100, mean = 7.5, sd = 1.8)
> > >prem_G2 <- rnorm(100, mean = 5, sd = 1.2)
> > >postsd_G2 <- rnorm(100, mean = 0.9, sd = 0.2)
> > >presd_G2 <- rnorm(100, mean = 0.6, sd = 0.2)
> > >n_G2 <- rpois(100, lambda = 50)
> > >
> > >G <- data.frame(postm_G1, prem_G1, postsd_G1, n_G1, presd_G1,
> > >postm_G2, prem_G2, postsd_G2, n_G2, presd_G2)
> > >
> > >###
> > >
> > >Thank you in advance for your time and help.
> > >
> > >Best wishes,
> > >Mika
> >
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-meta-analysis mailing list @ R-sig-meta-analysis using r-project.org
> To manage your subscription to this mailing list, go to:
> https://stat.ethz.ch/mailman/listinfo/r-sig-meta-analysis
>

	[[alternative HTML version deleted]]



More information about the R-sig-meta-analysis mailing list