[R-meta] Multilevel with dataset where lines are treatment means
Martineau, Roger
Roger.Martineau at AGR.GC.CA
Thu Sep 14 23:07:38 CEST 2017
Hi all,
I have a large data set to analyze and I want to get my model (and modus operandi) right before publication. Therefore I will provide a fictive example much smaller (n=28) but having a similar structure (i.e., experiment/study).
Thanks in advance,
Roger ☺
The fictive data are:
dat <- data.frame(
line=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28),
experiment=c("A","A","A","A","B","B","B","B","C","C","C","C",
"D","D","D","D","D","D","D","D","E","E","E","E","E","E","E","E"),
study=c("A1","A1","A2","A2","B1","B1","B2","B2","C1","C1","C1","C1",
"D1","D1","D1","D1","D2","D2","D2","D2","E1","E1","E1","E1","E2","E2","E2","E2"),
outcome=c(10,18,5,25,13,48,7,25,6,13,20,35,65,45,30,15,55,45,30,10,45,40,35,30,35,26,18,8),
SEM=c(5,5,6,6,7,7,3,3,7,7,7,7,6,6,6,6,5,5,5,5,6,6,6,6,5,5,5,5),
moderator1=c(5,10,25,45,8,30,5,35,5,10,15,20,10,20,25,30,20,30,40,50,10,15,20,25,5,6,9,12),
moderator2=c("cat1","cat1","cat1","cat1","cat1","cat1","cat1","cat1",
"cat1","cat1","cat1","cat1","cat2","cat2","cat2","cat2",
"cat2","cat2","cat2","cat2","cat2","cat2","cat2","cat2",
"cat2","cat2","cat2","cat2"))
dat
str(dat)
In this data set, "experiment" represents different publications (n=5) and one can fin 1 or 2 studies in each publication. Each study reports at least 2 treatment means (i.e., outcome) and up to 4 treatment means per study (see E2). In animal science, milk yield is an example of such an oucome. Indeed SEM of outcome is identical per study.
I included a continuous (moderator1) and a categorical moderator (moderator2). I deliberately set data such as there is an interaction (the outcome increase with increasing values of moderator1 for cat1 and vice-versa for cat2).
I have a couple of questions on the proper way to analyse such a data set reporting treatment means (and associated SEM).
**** Question 1 concerns the correct random argument structure to use:
# centering moderator:
dat <- within(dat, {
cmoderator1 <- moderator1 - mean(moderator1,na.rm = TRUE)
})
# model with 2 options for random argument:
res <- rma.mv(outcome,SEM^2, data=dat,
mods = ~
cmoderator1*factor(moderator2),
random = ~1|experiment/study, # option 1
#random = ~1|experiment/study/line, # option 2
method = "REML")
summary(res)
Results with option 1 and 2 are close but which one is correct ? and why ?
Based on previous posts, I assume that SE and P-values are correctly estimated using the robust function and the upper level cluster (i.e. experiment):
robust(res, cluster=dat$experiment)
In my large data set, I look for outliers and influential cases with studentitized residuals and Cook's distance values.
**** My second question is: is it OK to identify (and remove from the data set) outliers and influential cases on a study basis (i.e., using the cluster argument in the functions) ? or is it better to identify idividual treatment means (by removing the cluster argument) ?
#### .. Student residuals ####
tmp.stud.res <- rstudent(res,
cluster=dat$study, # is this OK ?
progbar=TRUE)
order(-tmp.stud.res$cluster$X2)
tmp.stud.res$cluster[c(xxx)]
#### .. Cook's distance ####
par(mfrow=c(1,1))
tmp.cook <- cooks.distance(res,
cluster=dat$study, # is this OK ?
progbar=TRUE)
plot(tmp.cook, type="o", pch=19)
which(tmp.cook > xxx)
Indeed, I do the other functions:
### Likelihood profile plots ###
par(mfrow=c(2,1))
profile.rma.mv(res)
### Normal Q-Q Plot
par(mfrow=c(1,1))
qqnorm(rstandard(res)$z, pch=19)
qqline(rstandard(res)$z)
### Funnel plots (publication bias) ###
par(mfrow=c(2,2))
funnel(res, main="Standard Error")
funnel(res, yaxis="vi", main="Sampling Variance")
funnel(res, yaxis="seinv", main="Inverse Standard Error")
funnel(res, yaxis="vinv", main="Inverse Sampling Variance")
This is the graph which shows the interaction:
#### .. Graph ####
(res <- rma.mv(outcome, SEM^2, data=dat,
mods = ~
factor(moderator2)-1 +
factor(moderator2):cmoderator1,
random = ~1|experiment/study,
method = "REML"))
tmp.coeff_cat <- round(coef(summary(res)), 4)
tmp <- ddply(dat, .(moderator2), summarise,
min = min(cmoderator1, na.rm = TRUE),
max = max(cmoderator1, na.rm = TRUE))
f1 <- function(x) { # function for tmp.casein studies
ifelse(x<tmp[1,3]& x>tmp[1,2],
tmp.coeff_cat[1,1] +
(tmp.coeff_cat[3,1]) *x,NA)}
f2 <- function(x) { # function for tmp.dietary studies
ifelse(x<tmp[2,3]& x>tmp[2,2],
tmp.coeff_cat[2,1] +
(tmp.coeff_cat[4,1]) *x,NA)}
tmp.plot <- ggplot(dat, aes(cmoderator1, outcome), na.rm = TRUE) +
geom_line(aes(group = study, colour = moderator2, linetype = moderator2), size = 1) +
geom_point(aes(group = study, colour = moderator2), size = 2)
tmp.plot + labs(
x = "Moderator1 (centered values)")+
scale_linetype_manual(
name="Category:",
values=c("solid","solid"),
breaks=c("cat1","cat2"),
labels=c("cat1","cat2")) +
scale_colour_manual(
name="Category:",
values=c("brown","royalblue"),
breaks=c("cat1","cat2"),
labels=c("cat1","cat2")) +
stat_function(fun=f1, colour = "red", size = 2, linetype="solid") +
stat_function(fun=f2, colour = "blue", size = 2, linetype="solid") +
ggtitle("Test")+
theme(
plot.title = element_text(hjust = 0.5, family="serif", size = 25),
axis.title=element_text(family= "serif", colour="black", size=18),
axis.line=element_line(color="black"),
axis.ticks = element_blank(),
axis.text=element_text(family= "serif", colour="black", size=15) ,
legend.title=element_text(family="serif", size=15),
legend.text=element_text(family="serif", size=15),
legend.position="bottom")
More information about the R-sig-meta-analysis
mailing list