The metafor Package

A Meta-Analysis Package for R

User Tools

Site Tools


tips:multiple_factors_interactions

Differences

This shows you the differences between two versions of the page.

Link to this comparison view

Both sides previous revisionPrevious revision
tips:multiple_factors_interactions [2022/08/03 11:35] Wolfgang Viechtbauertips:multiple_factors_interactions [2023/05/30 07:51] (current) Wolfgang Viechtbauer
Line 1: Line 1:
 ===== Models with Multiple Factors and Their Interaction ===== ===== Models with Multiple Factors and Their Interaction =====
  
-The example below shows how to test/examine multiple factors and their interaction in (mixed-effects) meta-regression models.+The example below shows how to test and examine multiple factors and their interaction in (mixed-effects) meta-regression models.
  
 ==== Data Preparation ==== ==== Data Preparation ====
  
 For the example, we will use the data from the meta-analysis by Raudenbush (1984) (see also Raudenbush & Bryk, 1985) of studies examining teacher expectancy effects on pupil IQ (''help(dat.raudenbush1985)'' provides a bit more background). The data can be loaded with: For the example, we will use the data from the meta-analysis by Raudenbush (1984) (see also Raudenbush & Bryk, 1985) of studies examining teacher expectancy effects on pupil IQ (''help(dat.raudenbush1985)'' provides a bit more background). The data can be loaded with:
 +
 <code rsplus> <code rsplus>
 library(metafor) library(metafor)
 dat <- dat.raudenbush1985 dat <- dat.raudenbush1985
 </code> </code>
 +
 I copy the dataset into ''dat'', which is a bit shorter and therefore easier to type further below. I copy the dataset into ''dat'', which is a bit shorter and therefore easier to type further below.
  
 For illustration purposes, we will categorize the ''weeks'' variable (i.e., the number of weeks of contact prior to the expectancy induction) into three levels (0 weeks = "none", 1-9 weeks = "some", and 10+ weeks = "high"): For illustration purposes, we will categorize the ''weeks'' variable (i.e., the number of weeks of contact prior to the expectancy induction) into three levels (0 weeks = "none", 1-9 weeks = "some", and 10+ weeks = "high"):
 +
 <code rsplus> <code rsplus>
 dat$weeks <- cut(dat$weeks, breaks=c(0,1,10,100), labels=c("none","some","high"), right=FALSE) dat$weeks <- cut(dat$weeks, breaks=c(0,1,10,100), labels=c("none","some","high"), right=FALSE)
Line 39: Line 42:
 19    19             Ginsburg 1970  some   group  aware -0.07 0.0303 19    19             Ginsburg 1970  some   group  aware -0.07 0.0303
 </code> </code>
 +
 Variable ''yi'' provides the standardized mean differences and ''vi'' the corresponding sampling variances for these 19 studies. For this illustration, we will focus on variables ''weeks'' (as described above) and ''tester'' (whether the test administrator was aware or blind to the researcher-provided designations of the children's intellectual potential) as possible moderator variables. Variable ''yi'' provides the standardized mean differences and ''vi'' the corresponding sampling variances for these 19 studies. For this illustration, we will focus on variables ''weeks'' (as described above) and ''tester'' (whether the test administrator was aware or blind to the researcher-provided designations of the children's intellectual potential) as possible moderator variables.
  
 The ''cut()'' command as used above automatically turns the ''weeks'' variable directly into a factor. Moreover, the ''none'' level is automatically set as the reference level (followed by levels ''some'' and ''high''). To turn the ''tester'' variable into a factor and to ensure that the ''blind'' level becomes the reference level in the models to be fitted below, we can use the following command: The ''cut()'' command as used above automatically turns the ''weeks'' variable directly into a factor. Moreover, the ''none'' level is automatically set as the reference level (followed by levels ''some'' and ''high''). To turn the ''tester'' variable into a factor and to ensure that the ''blind'' level becomes the reference level in the models to be fitted below, we can use the following command:
 +
 <code rsplus> <code rsplus>
 dat$tester <- relevel(factor(dat$tester), ref="blind") dat$tester <- relevel(factor(dat$tester), ref="blind")
Line 47: Line 52:
  
 A table with the number of factor level combinations (and margin counts) can be obtained with: A table with the number of factor level combinations (and margin counts) can be obtained with:
 +
 <code rsplus> <code rsplus>
 addmargins(table(dat$weeks, dat$tester)) addmargins(table(dat$weeks, dat$tester))
Line 57: Line 63:
   Sum      9    10  19   Sum      9    10  19
 </code> </code>
 +
 So, all factor level combinations are present in this dataset. So, all factor level combinations are present in this dataset.
  
Line 62: Line 69:
  
 We start by fitting a model that only contains the main effects of the two factors: We start by fitting a model that only contains the main effects of the two factors:
 +
 <code rsplus> <code rsplus>
 res.a1 <- rma(yi, vi, mods = ~ weeks + tester, data=dat) res.a1 <- rma(yi, vi, mods = ~ weeks + tester, data=dat)
Line 78: Line 86:
 QE(df = 15) = 24.1362, p-val = 0.0628 QE(df = 15) = 24.1362, p-val = 0.0628
  
-Test of Moderators (coefficient(s) 2,3,4):+Test of Moderators (coefficients 2:4):
 QM(df = 3) = 9.9698, p-val = 0.0188 QM(df = 3) = 9.9698, p-val = 0.0188
  
Line 90: Line 98:
  
 --- ---
-Signif. codes: ***’ 0.001 **’ 0.01 *’ 0.05 .’ 0.1 ‘ ’ 1+Signif. codes: '***0.001 '**0.01 '*0.05 '.0.1 ' ' 1
 </code> </code>
  
-Let's start with the table at the bottom containing the model results. The first coefficient in the table is the model intercept ($b_0 = .4020$), which is the estimated average standardized mean difference for levels ''none'' and ''blind'' of the ''weeks'' and ''tester'' factors, respectively. For this combination of levels, the effect is significantly different from zero ($p = .0020$), although the 95% confidence interval ($.1472$ to $.6567$) indicates quite some uncertainty how large the effect really is.+Let's start with the table at the bottom containing the model results. The first coefficient in the table is the model intercept ($b_0 = 0.4020$), which is the estimated average standardized mean difference for levels ''none'' and ''blind'' of the ''weeks'' and ''tester'' factors, respectively. For this combination of levels, the effect is significantly different from zero ($p = .0020$), although the 95% confidence interval ($0.1472$ to $0.6567$) indicates quite some uncertainty how large the effect really is.
  
-The second and third coefficients correspond to the dummy variables ''weekssome'' and ''weekshigh'' ($b_1 = -.2893$ and $b_2 = -.4422$, respectively). These coefficients estimate how much higher/lower the average effect is for these two levels of the ''weeks'' factorcompared to the reference level (i.e., ''none''). Therefore, these coefficients directly indicate the contrast between the ''none'' and the other two levels. Both coefficients are significantly different from zero ($p = .0334$ and $p = .0025$, respectively), so the average effect is significantly different (and in fact lower) when teachers had some or a high amount of contact with the students prior to the expectancy induction (instead of no contact).+The second and third coefficients correspond to the dummy variables ''weekssome'' and ''weekshigh'' ($b_1 = -0.2893$ and $b_2 = -0.4422$, respectively). These coefficients estimate how much higher/lower the average effect is for these two levels of the ''weeks'' factor compared to the reference level (i.e., ''none''). Therefore, these coefficients directly indicate the contrast between the ''none'' and the other two levels. Both coefficients are significantly different from zero ($p = .0334$ and $p = .0025$, respectively), so the average effect is significantly different (and in fact lower) when teachers had some or a high amount of contact with the students prior to the expectancy induction (instead of no contact).
  
 Note that we have not yet tested the ''weeks'' factor as a whole. To do so, we have to test the null hypothesis that both coefficients are equal to zero simultaneously (see also the illustration of how to [[tips:testing_factors_lincoms|test factors and linear combinations of parameters]] in meta-regression models). We can do this with: Note that we have not yet tested the ''weeks'' factor as a whole. To do so, we have to test the null hypothesis that both coefficients are equal to zero simultaneously (see also the illustration of how to [[tips:testing_factors_lincoms|test factors and linear combinations of parameters]] in meta-regression models). We can do this with:
 +
 <code rsplus> <code rsplus>
 anova(res.a1, btt=2:3) anova(res.a1, btt=2:3)
Line 105: Line 114:
 QM(df = 2) = 9.1550, p-val = 0.0103 QM(df = 2) = 9.1550, p-val = 0.0103
 </code> </code>
 +
 This is a Wald-type chi-square test with 2 degrees of freedom, indicating that the factor is in fact significant ($p = .0103$). This is a Wald-type chi-square test with 2 degrees of freedom, indicating that the factor is in fact significant ($p = .0103$).
  
 We have also not yet tested whether there is a difference between levels ''some'' and ''high'' of this factor. One could change the reference level of the ''weeks'' factor and refit the model to obtain this test. Alternatively, and more elegantly, we can just test the difference between these two coefficients directly. We can do this with: We have also not yet tested whether there is a difference between levels ''some'' and ''high'' of this factor. One could change the reference level of the ''weeks'' factor and refit the model to obtain this test. Alternatively, and more elegantly, we can just test the difference between these two coefficients directly. We can do this with:
 +
 <code rsplus> <code rsplus>
 anova(res.a1, X=c(0,1,-1,0)) anova(res.a1, X=c(0,1,-1,0))
Line 124: Line 135:
  
 The ''linearHypothesis()'' function from the ''car'' package also allows us to conduct such a test. After (installing and) loading the package, either one of the following commands leads to the same results: The ''linearHypothesis()'' function from the ''car'' package also allows us to conduct such a test. After (installing and) loading the package, either one of the following commands leads to the same results:
 +
 <code rsplus> <code rsplus>
 library(car) library(car)
Line 142: Line 154:
 2     15  1 2.265     0.1323 2     15  1 2.265     0.1323
 </code> </code>
-So, we actually do not have enough evidence to conclude that there is a significant difference between these two levels ($p = .1323$).+ 
 +So, we actually do not have sufficient evidence to conclude that there is a significant difference between these two levels ($p = .1323$).
  
 The ''glht()'' function from the ''multcomp'' package also allows for such tests and actually makes it easy to conduct all pairwise comparisons between factor levels (with or without adjusted p-values due to multiple testing). To test all three linear combinations against each other, we would use: The ''glht()'' function from the ''multcomp'' package also allows for such tests and actually makes it easy to conduct all pairwise comparisons between factor levels (with or without adjusted p-values due to multiple testing). To test all three linear combinations against each other, we would use:
 +
 <code rsplus> <code rsplus>
 library(multcomp) library(multcomp)
Line 161: Line 175:
 (Adjusted p values reported -- none method) (Adjusted p values reported -- none method)
 </code> </code>
 +
 No adjustment to the p-values has been used for these tests (''test=adjusted("none")''), but other options are possible for the ''test'' argument (see ''help(summary.glht)'' and the illustration of how to [[tips:testing_factors_lincoms|test factors and linear combinations of parameters]]). Note that the p-values for the three tests above are exactly the same as the corresponding two p-values from the model table (coefficients 2 and 3) and the one obtained with the ''linearHypothesis()'' function above. No adjustment to the p-values has been used for these tests (''test=adjusted("none")''), but other options are possible for the ''test'' argument (see ''help(summary.glht)'' and the illustration of how to [[tips:testing_factors_lincoms|test factors and linear combinations of parameters]]). Note that the p-values for the three tests above are exactly the same as the corresponding two p-values from the model table (coefficients 2 and 3) and the one obtained with the ''linearHypothesis()'' function above.
  
-Finally, the last (4th) coefficient in the model table corresponds to the dummy variable ''testeraware'' ($b_3 = -.0511$) and estimates how much higher/lower the average effect is for studies where the test administrator was aware of the children's designationcompared to studies where the test administrator was blind. So, again, this is a contrast between these two levels. The coefficient is not significantly different from zero ($p = .5815$), so this factor does not appear to account for (some of the) heterogeneity in the effects.+Finally, the last (4th) coefficient in the model table corresponds to the dummy variable ''testeraware'' ($b_3 = -0.0511$) and estimates how much higher/lower the average effect is for studies where the test administrator was aware of the children's designation compared to studies where the test administrator was blind. So, again, this is a contrast between these two levels. The coefficient is not significantly different from zero ($p = .5815$), so this factor does not appear to account for (some of the) heterogeneity in the effects.
  
 The output includes some additional results. First of all, the results given under ''Test of Moderators (coefficient(s) 2,3,4)'' is the [[wp>Omnibus_test|omnibus test]] of coefficients 2, 3, and 4 (the first being the intercept, which is not included in the test by default). This Wald-type chi-square test is significant ($Q_M = 9.97, df = 3, p = .02$), so we can reject the null hypothesis that all of these coefficients are equal to zero (which would imply that there is no relationship between the size of the effect and these two factors). Moreover, the pseudo-$R^2$ value given suggests that roughly 50% of the total amount of heterogeneity is accounted for when including these two factors in the model. Finally, the ''Test for Residual Heterogeneity'' just falls short of being significant ($Q_E = 24.14, df = 15, p = .06$), which would suggest that there is no significant amount of residual heterogeneity in the true effects when including these two factors in the model. The output includes some additional results. First of all, the results given under ''Test of Moderators (coefficient(s) 2,3,4)'' is the [[wp>Omnibus_test|omnibus test]] of coefficients 2, 3, and 4 (the first being the intercept, which is not included in the test by default). This Wald-type chi-square test is significant ($Q_M = 9.97, df = 3, p = .02$), so we can reject the null hypothesis that all of these coefficients are equal to zero (which would imply that there is no relationship between the size of the effect and these two factors). Moreover, the pseudo-$R^2$ value given suggests that roughly 50% of the total amount of heterogeneity is accounted for when including these two factors in the model. Finally, the ''Test for Residual Heterogeneity'' just falls short of being significant ($Q_E = 24.14, df = 15, p = .06$), which would suggest that there is no significant amount of residual heterogeneity in the true effects when including these two factors in the model.
  
 Based on the fitted model, we can also obtain the estimated average effect for each combination of the factor levels. For this, we can use the ''predict()'' function, specifying how the dummy variables should be filled in. For levels ''none'', ''some'', and ''high'' of the ''weeks'' factor and level ''blind'' of the ''tester'' factor, we would use: Based on the fitted model, we can also obtain the estimated average effect for each combination of the factor levels. For this, we can use the ''predict()'' function, specifying how the dummy variables should be filled in. For levels ''none'', ''some'', and ''high'' of the ''weeks'' factor and level ''blind'' of the ''tester'' factor, we would use:
 +
 <code rsplus> <code rsplus>
 predict(res.a1, newmods=rbind(c(0,0,0),c(1,0,0),c(0,1,0))) predict(res.a1, newmods=rbind(c(0,0,0),c(1,0,0),c(0,1,0)))
Line 177: Line 193:
 3 -0.0402 0.1020 -0.2401 0.1597 -0.3163 0.2359 3 -0.0402 0.1020 -0.2401 0.1597 -0.3163 0.2359
 </code> </code>
 +
 For the same levels of the ''weeks'' factor and level ''aware'' of the ''tester'' factor, we would use: For the same levels of the ''weeks'' factor and level ''aware'' of the ''tester'' factor, we would use:
 +
 <code rsplus> <code rsplus>
 predict(res.a1, newmods=rbind(c(0,0,1),c(1,0,1),c(0,1,1))) predict(res.a1, newmods=rbind(c(0,0,1),c(1,0,1),c(0,1,1)))
Line 187: Line 205:
 3 -0.0913 0.0857 -0.2592 0.0767 -0.3452 0.1627 3 -0.0913 0.0857 -0.2592 0.0767 -0.3452 0.1627
 </code> </code>
 +
 Note that, by default, the intercept is automatically included in the calculation of these predicted values (so only a vector of length 3 or a matrix with 3 columns should be specified via the ''newmods'' argument). The values under ''ci.lb'' and ''ci.ub'' are the bounds of the 95% confidence intervals, while the values under ''pi.lb'' and ''pi.ub'' are the bounds of the 95% prediction intervals (see ''help(predict.rma)'' for more details). Note that, by default, the intercept is automatically included in the calculation of these predicted values (so only a vector of length 3 or a matrix with 3 columns should be specified via the ''newmods'' argument). The values under ''ci.lb'' and ''ci.ub'' are the bounds of the 95% confidence intervals, while the values under ''pi.lb'' and ''pi.ub'' are the bounds of the 95% prediction intervals (see ''help(predict.rma)'' for more details).
  
 We can use an alternative model specification, where we leave out the model intercept: We can use an alternative model specification, where we leave out the model intercept:
 +
 <code rsplus> <code rsplus>
 res.a2 <- rma(yi, vi, mods = ~ weeks + tester - 1, data=dat) res.a2 <- rma(yi, vi, mods = ~ weeks + tester - 1, data=dat)
Line 219: Line 239:
 Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 </code> </code>
 +
 Note that the first three coefficients now directly correspond to the estimated average effects for levels ''none'', ''some'', and ''high'' of the ''weeks'' factor and level ''blind'' of the ''tester'' factor. The last (4th) coefficient again indicates how much higher/lower the average effect is when the level of the ''tester'' factor is changed to ''aware''. Also, note that the omnibus test now includes all four coefficients by default and has a different interpretation: It is essentially testing the null hypothesis that the true average effect is zero across all factor level combinations. Note that the first three coefficients now directly correspond to the estimated average effects for levels ''none'', ''some'', and ''high'' of the ''weeks'' factor and level ''blind'' of the ''tester'' factor. The last (4th) coefficient again indicates how much higher/lower the average effect is when the level of the ''tester'' factor is changed to ''aware''. Also, note that the omnibus test now includes all four coefficients by default and has a different interpretation: It is essentially testing the null hypothesis that the true average effect is zero across all factor level combinations.
  
 The model including only the main effects of the two factors assumes that the influence of the two factors is additive. We can illustrate the implications graphically by plotting the estimated average effects for all factor level combinations: The model including only the main effects of the two factors assumes that the influence of the two factors is additive. We can illustrate the implications graphically by plotting the estimated average effects for all factor level combinations:
 +
 <code rsplus> <code rsplus>
 plot(coef(res.a2)[1:3], type="o", pch=19, xlim=c(0.8,3.2), ylim=c(-0.1,0.5), xlab="Amount of Contact Prior to Expectancy Induction", ylab="Standardized Mean Difference", xaxt="n", bty="l") plot(coef(res.a2)[1:3], type="o", pch=19, xlim=c(0.8,3.2), ylim=c(-0.1,0.5), xlab="Amount of Contact Prior to Expectancy Induction", ylab="Standardized Mean Difference", xaxt="n", bty="l")
Line 237: Line 259:
  
 Now let's move on to a model where we allow the two factors in the model to interact (note that it is questionable whether such a complex model should even be considered for these data; again, these analyses are for illustration purposes only). Such a model can be fitted with: Now let's move on to a model where we allow the two factors in the model to interact (note that it is questionable whether such a complex model should even be considered for these data; again, these analyses are for illustration purposes only). Such a model can be fitted with:
 +
 <code rsplus> <code rsplus>
 res.i1 <- rma(yi, vi, mods = ~ weeks * tester, data=dat) res.i1 <- rma(yi, vi, mods = ~ weeks * tester, data=dat)
Line 270: Line 293:
 </code> </code>
  
-As before, the first coefficient in the table is the model intercept ($b_0 = .3067$), which is the estimated average effect for levels ''none'' and ''blind'' of the ''weeks'' and ''tester'' factors, respectively. The second and third coefficients ($b_1 = -.1566$ and $b_2 = -.3267$, respectively) indicate how much higher/lower the estimated average effect is for levels ''some'' and ''high'' of the ''weeks'' factor when the level of the ''tester'' factor is held at ''blind''. The fourth coefficient ($b_3 = .1759$) provides an estimate of how much higher/lower the average effect is for the ''aware'' level of the ''tester'' factor when the level of the ''weeks'' factor is held at ''none''. Finally, the last two coefficients ($b_4 = -.2775$ and $b_5 = -.2634$, respectively) indicate how much higher/lower the contrast between the ''aware'' and ''blind'' levels of the ''tester'' factor are for levels ''some'' and ''high'' of the ''weeks'' factor (or equivalently, how much higher/lower the contrasts between the ''some'' and the ''none'' levels and the ''high'' and the ''none'' levels of the ''weeks'' factor are for level ''aware'' of the ''tester'' factor).+As before, the first coefficient in the table is the model intercept ($b_0 = 0.3067$), which is the estimated average effect for levels ''none'' and ''blind'' of the ''weeks'' and ''tester'' factors, respectively. The second and third coefficients ($b_1 = -0.1566$ and $b_2 = -0.3267$, respectively) indicate how much higher/lower the estimated average effect is for levels ''some'' and ''high'' of the ''weeks'' factor when the level of the ''tester'' factor is held at ''blind''. The fourth coefficient ($b_3 = 0.1759$) provides an estimate of how much higher/lower the average effect is for the ''aware'' level of the ''tester'' factor when the level of the ''weeks'' factor is held at ''none''. Finally, the last two coefficients ($b_4 = -0.2775$ and $b_5 = -0.2634$, respectively) indicate how much higher/lower the contrast between the ''aware'' and ''blind'' levels of the ''tester'' factor are for levels ''some'' and ''high'' of the ''weeks'' factor (or equivalently, how much higher/lower the contrasts between the ''some'' and the ''none'' levels and the ''high'' and the ''none'' levels of the ''weeks'' factor are for level ''aware'' of the ''tester'' factor).
  
 These last two coefficients allow the influence of the ''aware'' factor to differ across the various levels of the ''weeks'' factor (and vice-versa). Therefore, if these two coefficients are in reality both equal to zero, then there is in fact no interaction. So, in order to test whether an interaction is actually present, we just need to test these two coefficients, which we can do with: These last two coefficients allow the influence of the ''aware'' factor to differ across the various levels of the ''weeks'' factor (and vice-versa). Therefore, if these two coefficients are in reality both equal to zero, then there is in fact no interaction. So, in order to test whether an interaction is actually present, we just need to test these two coefficients, which we can do with:
 +
 <code rsplus> <code rsplus>
 anova(res.i1, btt=5:6) anova(res.i1, btt=5:6)
Line 280: Line 304:
 QM(df = 2) = 0.8576, p-val = 0.6513 QM(df = 2) = 0.8576, p-val = 0.6513
 </code> </code>
 +
 Based on this Wald-type chi-square test, we find no significant evidence that an interaction is actually present ($Q_M = 0.86, df = 2, p = .65$). Based on this Wald-type chi-square test, we find no significant evidence that an interaction is actually present ($Q_M = 0.86, df = 2, p = .65$).
  
 Alternatively, we could use a likelihood ratio test (LRT) to test the interaction. For this, we must fit both models with maximum likelihood (ML) estimation and then conduct a model comparison: Alternatively, we could use a likelihood ratio test (LRT) to test the interaction. For this, we must fit both models with maximum likelihood (ML) estimation and then conduct a model comparison:
 +
 <code rsplus> <code rsplus>
 res.a1.ml <- rma(yi, vi, mods = ~ weeks + tester, data=dat, method="ML") res.a1.ml <- rma(yi, vi, mods = ~ weeks + tester, data=dat, method="ML")
Line 289: Line 315:
 </code> </code>
 <code output> <code output>
-        df    AIC     BIC    AICc logLik    LRT   pval      QE  tau^2    R^2+        df    AIC     BIC    AICc logLik    LRT   pval      QE  tau^2      R^2
 Full     7 8.5528 15.1639 18.7346 2.7236               23.4646 0.0000 Full     7 8.5528 15.1639 18.7346 2.7236               23.4646 0.0000
-Reduced  5 5.2247  9.9469  9.8401 2.3876 0.6719 0.7147 24.1362 0.0000 71.65%+Reduced  5 5.2247  9.9469  9.8401 2.3876 0.6719 0.7147 24.1362 0.0000 71.6487%
 </code> </code>
 +
 The results from the LRT ($\chi^2 = 0.67, df = 2, p = .71$) lead to the same conclusion. The results from the LRT ($\chi^2 = 0.67, df = 2, p = .71$) lead to the same conclusion.
  
 For illustration purposes, let us however proceed with a further examination of this model. As before, the ''anova()'', ''linearHypothesis()'', and ''glht()'' functions could be used to test factor levels against each other and to obtain pairwise comparisons. Also, the ''predict()'' function can be used to obtain the estimated average effect for particular factor level combinations. For example, to obtain the estimated average effect for level ''high'' of the ''weeks'' factor and level ''aware'' of the ''tester'' factor, we would use: For illustration purposes, let us however proceed with a further examination of this model. As before, the ''anova()'', ''linearHypothesis()'', and ''glht()'' functions could be used to test factor levels against each other and to obtain pairwise comparisons. Also, the ''predict()'' function can be used to obtain the estimated average effect for particular factor level combinations. For example, to obtain the estimated average effect for level ''high'' of the ''weeks'' factor and level ''aware'' of the ''tester'' factor, we would use:
 +
 <code rsplus> <code rsplus>
 predict(res.i1, newmods=c(0,1,1,0,1)) predict(res.i1, newmods=c(0,1,1,0,1))
Line 305: Line 333:
  
 We can also make use of an alternative parameterization of the model, which is particularly useful in this context: We can also make use of an alternative parameterization of the model, which is particularly useful in this context:
 +
 <code rsplus> <code rsplus>
 res.i2 <- rma(yi, vi, mods = ~ weeks:tester - 1, data=dat) res.i2 <- rma(yi, vi, mods = ~ weeks:tester - 1, data=dat)
Line 338: Line 367:
  
 Now, the table with the model results directly provides the estimated average effect for each factor level combination. It is now also quite easy to test particular factor level combinations against each other. For example, to test the difference between levels ''some'' and ''high'' of the ''weeks'' factor within the ''blind'' level of the ''tester'' factor, we could use: Now, the table with the model results directly provides the estimated average effect for each factor level combination. It is now also quite easy to test particular factor level combinations against each other. For example, to test the difference between levels ''some'' and ''high'' of the ''weeks'' factor within the ''blind'' level of the ''tester'' factor, we could use:
 +
 <code rsplus> <code rsplus>
 anova(res.i2, X=c(0,1,-1,0,0,0)) anova(res.i2, X=c(0,1,-1,0,0,0))
Line 354: Line 384:
  
 Or with the ''linearHypothesis()'' function (abbreviated output shown): Or with the ''linearHypothesis()'' function (abbreviated output shown):
 +
 <code rsplus> <code rsplus>
 linearHypothesis(res.i2, c("weekssome:testerblind - weekshigh:testerblind = 0")) linearHypothesis(res.i2, c("weekssome:testerblind - weekshigh:testerblind = 0"))
Line 364: Line 395:
  
 To test the same contrast within the ''aware'' level of the ''tester'' factor, we would use: To test the same contrast within the ''aware'' level of the ''tester'' factor, we would use:
 +
 <code rsplus> <code rsplus>
 anova(res.i2, X=c(0,0,0,0,1,-1)) anova(res.i2, X=c(0,0,0,0,1,-1))
Line 380: Line 412:
  
 Or with the ''linearHypothesis()'' function (abbreviated output shown): Or with the ''linearHypothesis()'' function (abbreviated output shown):
 +
 <code rsplus> <code rsplus>
 linearHypothesis(res.i2, c("weekssome:testeraware - weekshigh:testeraware = 0")) linearHypothesis(res.i2, c("weekssome:testeraware - weekshigh:testeraware = 0"))
Line 389: Line 422:
 </code> </code>
  
-If one were inclined to do so, one could even test all 6 factor level combinations against each other (i.e., all possible pairwise comparisons). The ''glht()'' function combined with the ''contMat()'' function makes this relatively easy:+If one were inclined to do so, one could even test all 6 factor level combinations against each other (i.e., all possible pairwise comparisons). The ''glht()'' function combined with the ''contrMat()'' function makes this relatively easy: 
 <code rsplus> <code rsplus>
 summary(glht(res.i2, linfct=contrMat(c("none:blind"=1,"some:blind"=1,"high:blind"=1, summary(glht(res.i2, linfct=contrMat(c("none:blind"=1,"some:blind"=1,"high:blind"=1,
Line 420: Line 454:
 (Adjusted p values reported -- none method) (Adjusted p values reported -- none method)
 </code> </code>
 +
 As before, no adjustment to the p-values has been used for these tests, but this would probably be advisable when conducting so many tests. As before, no adjustment to the p-values has been used for these tests, but this would probably be advisable when conducting so many tests.
  
 We can also illustrate the results from the interaction model by plotting the estimated average effects for all factor level combinations: We can also illustrate the results from the interaction model by plotting the estimated average effects for all factor level combinations:
 +
 <code rsplus> <code rsplus>
 plot(coef(res.i2)[1:3], type="o", pch=19, xlim=c(0.8,3.2), ylim=c(-0.1,0.5), xlab="Amount of Contact Prior to Expectancy Induction", ylab="Standardized Mean Difference", xaxt="n", bty="l") plot(coef(res.i2)[1:3], type="o", pch=19, xlim=c(0.8,3.2), ylim=c(-0.1,0.5), xlab="Amount of Contact Prior to Expectancy Induction", ylab="Standardized Mean Difference", xaxt="n", bty="l")
tips/multiple_factors_interactions.txt · Last modified: 2023/05/30 07:51 by Wolfgang Viechtbauer