[R] matafor package - categorical moderator interpretation question
Viechtbauer Wolfgang (SP)
wolfgang.viechtbauer at maastrichtuniversity.nl
Tue Apr 4 08:41:28 CEST 2017
You are not conducting a proper test of the moderator. When you use 'mods = ~factor(sample_data$Participants) - 1', the model does not include an intercept term but dummy variables corresponding to all levels of the moderator. The omnibus test you are getting therefore tests the null hypothesis that the model coefficients corresponding to the dummy variables are all simultaenously equal to zero. What you want to do is test the null hypothesis that the coefficients are equal to each other. The easiest way to obtain this test is to use 'mods = ~factor(sample_data$Participants)'.
Best,
Wolfgang
--
Wolfgang Viechtbauer, Ph.D., Statistician | Department of Psychiatry and
Neuropsychology | Maastricht University | P.O. Box 616 (VIJV1) | 6200 MD
Maastricht, The Netherlands | +31 (43) 388-4170 | http://www.wvbauer.com
-----Original Message-----
From: R-help [mailto:r-help-bounces at r-project.org] On Behalf Of Calin-Jageman, Robert
Sent: Monday, April 03, 2017 23:32
To: r-help at r-project.org
Subject: [R] matafor package - categorical moderator interpretation question
What does it mean if a categorical moderator is significant overall but has no significant pairwise contrasts between moderator levels?
I'm using metaphor to conduct a meta-analysis with a categorical moderator with 3 levels; this yields a significant result:
Test of Moderators (coefficient(s) 1,2,3):
F(df1 = 3, df2 = 37) = 4.6052, p-val = 0.0078
Model Results:
estimate se tval pval ci.lb ci.ub
factor(sample_data$Participants)Adults 0.3920 0.2847 1.3771 0.1768 -0.1848 0.9688
factor(sample_data$Participants)Online 0.1403 0.1283 1.0935 0.2812 -0.1197 0.4004
factor(sample_data$Participants)Students 0.2350 0.0717 3.2747 0.0023 0.0896 0.3803 **
But then I conduct contrasts between each moderator level, and none of these are significant (no correction for multiple comparisons applied):
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
Online - Adults == 0 -0.2517 0.3123 -0.806 0.420
Students - Adults == 0 -0.1571 0.2936 -0.535 0.593
Students - Online == 0 0.0946 0.1470 0.643 0.520
(Adjusted p values reported -- none method)
Any thoughts or guides to interpretation are appreciated! My code and sample data are at the end of the email. My interpretation is that while one of the moderator levels may have be a significant factor in the overall analysis, the comparisons between moderator levels are noisier because they test to see if there is a difference in the weights between the two levels. Given this pattern of results, I conclude the different moderator levels are probably not strong predictors of effect size. I'm a bit uncertain if this is correct, and would appreciate any feedback.
Bob
========
Robert Calin-Jageman
Professor, Psychology
Neuroscience Program Director
Dominican University
Parmer 210
7900 West Division
River Forest, IL 60305
rcalinjageman at dom.edu
708.524.6581
http://calin-jageman.net
Sample data link:
https://www.dropbox.com/s/hzz9wmt1d9tcxsm/red_effect_males.csv?dl=0
Code:
#load required libraries
library("metafor")
library("multcomp")
sample_data <- read.csv("red_effect_males.csv")
#Overall test of categorical moderator, reports significant result
mod_test = rma(yi, vi, mods = ~factor(sample_data$Participants) - 1, data=sample_data, knha = TRUE)
print(mod_test)
#Now do pairwise contrasts - but these show no significant contrasts....why?
cont_holder <- c(1:length(unique(sample_data$Participants)))
names(cont_holder) <- sort(unique(sample_data$Participants))
print(summary(glht(mod_test, linfct=contrMat(cont_holder, "Tukey")), test=adjusted("none")))
#Now print individual meta-analysis for each subgroub... Effect sizes estimates and CIs aren't the same as in overall analysis...why?
subgroup_list <- split(sample_data, sample_data$Participants, drop=FALSE)
for (subgroup in subgroup_list) {
print(paste("Individual results for: ", subgroup$Participants[1]))
print(rma(yi, vi, data=subgroup, knha=TRUE))
}
More information about the R-help
mailing list