4

Consider three different studies based on three different datasets. Each of them has a continuous predictor x, a dichotomous predictor z, and a continuous outcome y.

In each of these three studies, researchers analyze the x*z interaction on y.

Here are three reproducible datasets that meet these criteria:

# generating data
set.seed(1839) # set seed
# dataset 1
dat1 <- data.frame(x=rnorm(50), z=factor(c(rep("A", 25), rep("B", 25))))
dat1$y <- c(rnorm(25), dat1$x[26:50]+rnorm(25,0,2))
# dataset 2
dat2 <- data.frame(x=rnorm(50), z=factor(c(rep("A", 25), rep("B", 25))))
dat2$y <- c(rnorm(25), dat2$x[26:50]+rnorm(25,0,2.3))
# dataset 3
dat3 <- data.frame(x=rnorm(50), z=factor(c(rep("A", 25), rep("B", 25))))
dat3$y <- c(rnorm(25), dat3$x[26:50]+rnorm(25,0,1.5))

Now, here are the coefficient tables from each of the three interactions:

> # data 1 results
> summary(lm(y~x*z, dat1))$coef
               Estimate Std. Error     t value   Pr(>|t|)
(Intercept) -0.27501165  0.2685408 -1.02409633 0.31114537
x            0.02228078  0.3083321  0.07226228 0.94270647
zB          -0.59286879  0.3791864 -1.56352861 0.12478258
x:zB         0.70988125  0.3978048  1.78449621 0.08093913

> # data 2 results
> summary(lm(y~x*z, dat2))$coef
               Estimate Std. Error    t value    Pr(>|t|)
(Intercept)  0.09181368  0.3377780  0.2718166 0.786979289
x           -0.31598416  0.3566014 -0.8860990 0.380173233
zB          -0.06912280  0.4773531 -0.1448043 0.885497964
x:zB         1.66773706  0.4983741  3.3463555 0.001638411

> # data 3 results
> summary(lm(y~x*z, dat3))$coef
               Estimate Std. Error     t value  Pr(>|t|)
(Intercept)  0.14508485  0.2977212  0.48731777 0.6283475
x            0.03496345  0.3963061  0.08822337 0.9300821
zB          -0.34413335  0.4124389 -0.83438616 0.4083758
x:zB         0.30622182  0.4672678  0.65534538 0.5155096

Is there a way to get a meta-analytic estimate of x:zB across these three studies?

Other things to consider:

  1. I have access to the raw data, and the scales x and y are all the same, but the manipulations involved in z were slightly different (although conceptually the same). I feel as if it is not best practice to simply collapse across the three datasets; is there any support for my intuition here? Or would collapsing across the datasets be defensible?

  2. I know that one could always get the effect of z as a Cohen's d, meta-analyze that, and use the mean of x from each study as predictors in a meta-regression. But note that the small number of studies here (three) makes that untenable.

Mark White
  • 10,252
  • 1
    Would you consider obtaining the separate correlations between x and y for each level of z? Then you could easily create a data frame that you could meta-analyze with a vector of correlations, a vector of their variances, and a vector of values indicator group (in your example A or B). – Matt Barstead Jul 12 '17 at 02:20
  • So a meta-analysis with six correlations, and use group (A or B) as a predictor in meta-regression (or subgroup analysis)? – Mark White Jul 12 '17 at 02:25
  • That would be my instinct, depending on how tenable/untenable it is to just collapse across the different z's in each study. I think your intuition expressed in #1 is probably on point, though. At the very least it might be hard to defend to a reviewer. In an ideal world you would have a few more correlations, but when do we ever get everything we need from our data? I'd be curious to see what other responses you get in this case. – Matt Barstead Jul 12 '17 at 02:35
  • 2
    Since you have the raw data, how about you run it as a single model adding dummy coded variables to represent the study the data come from? You could inflate the standard errors to account for clustering using the square root of the design effect. This could be a fixed effect approach to accounting for clustering when the number of clusters is small. – Heteroskedastic Jim Jul 12 '17 at 02:42
  • @user162986 what about a multilevel model? Nest people within study: lmer(y ~ x*z + (1+x+z|studyid) – Mark White Jul 12 '17 at 13:42
  • 1
    The number of clusters you have would be too small for any of the estimators to work correctly. This has been studied by Dan McNeish and Laura Stapleton, and they recommend this fixed effects approach when the number of clusters is this small. Their simulation only went as low as 4. McNeish, D., & Stapleton, L. M. (2016). Modeling Clustered Data with Very Few Clusters. Multivariate Behavioral Research, 51(4), 495–518. https://doi.org/10.1080/00273171.2016.1167008 – Heteroskedastic Jim Jul 12 '17 at 13:58

2 Answers2

2

I would conduct this as a single regression analysis since you have the raw data. However, the clustering of observations within studies would violate statistical independence, hence possibly deflating the standard errors of our regression coefficients.

To account for the clustering, the standard response might be to perform a mixed-effects analysis (a multi-level model). However, the number of clusters is very small (n=3). When this is the case, the estimators for multi-level models are all faulty. The same goes for other most methods of accounting for clustering - cluster-robust standard errors, ...

Bayesian methods could suffice. However, a fixed-effects approach is recommended by McNeish and Stapleton (2016) who performed simulations with as few as four clusters. To account for clustering, they recommend inflating the standard errors of our regression coefficients using the square root of the design effect.

Some sample code to do this:

# Continuing from OP's example
dat <- rbind(dat1, dat2, dat3)
dat <- cbind(dat, g=c(rep("A", 50), rep("B", 50), rep("C", 50)))
str(dat)
'data.frame':   150 obs. of  4 variables:
 $ x: num  1.013 -0.685 0.349 -1.625 -0.516 ...
 $ z: Factor w/ 2 levels "A","B": 1 1 1 1 1 1 1 1 1 1 ...
 $ y: num  0.686 -0.827 -0.507 0.117 0.504 ...
 $ g: Factor w/ 3 levels "A","B","C": 1 1 1 1 1 1 1 1 1 1 ...

# Calculate ICC, design effect and the root of the design effect
if (!require(ICC)) { install.packages("ICC") }
icc.est <- ICCest(g, y, dat)
icc <- icc.est$ICC
k <- icc.est$k # Average size of the clusters, here would be 50
deff <- 1 + icc * (k - 1)
deft <- sqrt(deff)

# Conduct regression and inflate standard errors
(model.lm <- lm(y ~ x*z + g, data = dat))

Call:
lm(formula = y ~ x * z + g, data = dat)

Coefficients:
(Intercept)            x           zB           gB           gC         x:zB  
    -0.3980      -0.1081      -0.3285       0.6191       0.5997       0.8671

model.lm <- summary(model.lm)
model.lm.se <- as.numeric(model.lm$coefficients[, 2]) # Obtain standard errors
model.lm.se <- model.lm.se * deft # Multiply them by DEFT
t <- qt(.975, model.lm$df[2]) # Obtain coefficient of se
lower.bound <- model.lm$coefficients[, 1] - t * model.lm.se
upper.bound <- model.lm$coefficients[, 1] + t * model.lm.se
(final.results <- data.frame(
  estimate = as.numeric(model.lm$coefficients[, 1]),
  se = model.lm.se, lb = lower.bound, up = upper.bound
))
              estimate        se          lb        up
(Intercept) -0.3980202 0.3793703 -1.14787409 0.3518336
x           -0.1080681 0.3146023 -0.72990321 0.5137670
zB          -0.3285170 0.3803951 -1.08039639 0.4233623
gB           0.6190885 0.4635929 -0.29723767 1.5354147
gC           0.5997314 0.4633536 -0.31612175 1.5155846
x:zB         0.8671220 0.4047123  0.06717783 1.6670662

final.results now contains your modeling results. You can see the estimates for the different predictors and their 95% confidence intervals.

McNeish, D., & Stapleton, L. M. (2016). Modeling Clustered Data with Very Few Clusters. Multivariate Behavioral Research, 51(4), 495–518. https://doi.org/10.1080/00273171.2016.1167008

Heteroskedastic Jim
  • 5,167
  • This is great, thanks! Why is g included as a predictor? Because that's where we get the ICC from, and it needs to be in there to account for clustering? It seems like it makes it a predictor, which the paper says FEM cannot do. – Mark White Jul 12 '17 at 22:39
  • g is included in the final model not to calculate the ICC, that is a separate calculation. It is present to include the clusters as fixed effects, dummy variables:

    "With FEMs (a.k.a. dummy variable regression), cluster affiliation indicators (0/1 indicator variables, one for each cluster in the data) are included in the model as predictor variables with the goal being to account for the nested structure of the data without estimating the random effects, particularly when assumptions inherent with random effects are untenable or estimation may be computationally complex"

     – Heteroskedastic Jim Jul 12 '17 at 22:43
  • Ah, gotcha. I'm a big MLM fan, so I'm having a hard time wrapping my head around using cluster variables as fixed effect predictors. – Mark White Jul 12 '17 at 22:50
  • Fixed effects does seem to be the way to go here as otherwise you are estimating a variance from three observations which must be fairly unstable as an estimate.. – mdewey Jul 13 '17 at 11:10
1

If you have access to the raw data then it would seem best to use that. You would want to include a three level factor to distinguish study and probably also interact it with the variables you are really interested in in case there is heterogeneity there.

I definitely would not do your second option as this turns it from a study of the effect of $x$ as an individual level predictor to the effect of $x$ as an ecological predictor which is not the same at all.

mdewey
  • 17,806
  • I might hold off that interaction. Since OP is interested in the interaction, then OP would then create a three-way interaction, and would have to probe it to assess its significance (not just statistical significance). I might be more comfortable with establishing that the qualitative conclusions reached from the interactions in the three studies are similar enough, such that I can assume there is no such heterogeneity. Then I would use the square root of the design effect to inflate the standard errors to account for clustered nature of the analysis. – Heteroskedastic Jim Jul 12 '17 at 14:10
  • @user162986 from the paper you cited, it seems like MLMs have problems with convergence (as a side note, I'm not sure how that speaks statistically to the issue, other than low convergence rates means we cannot use the model). However, I have no Level 2 predictors. The simulation study includes a Level 2 predictor and two cross-level interactions. In lme4::lmer syntax, the model would be: y~x*z + (1|g), where g is the cluster variable. Conceptually, why would an MLM still do a poor job with a design that simple? – Mark White Jul 12 '17 at 22:48
  • I think the sample size is problematic. It's like doing statistics with a sample size of 3. The coverage rates deviates from the nominal rate for MLM. I'm re-reading the paper and I am not sure I implemented the solution correctly. They seem to be concerned with level-2 predictors and may have recommended the correction for level-2 predictors only. My code may have been much ado about nothing. But there is nothing wrong using clusters as FE. It just limits the inference to the clusters in the study, as in a fixed-effects meta-analysis. – Heteroskedastic Jim Jul 12 '17 at 23:14