# Meta-analysis analysis script DN February 1st 2018, edited by GP on 3rd May 2018 # Data analysis script for telomere meta-analysis # Libraries library(metafor) library(ggplot2) library(psych) # Requires other packages for supplementary analyses, see relevant section # Read in data d=read.csv("processed data.csv") #Make reduced dataset with just the sample sizes greater than 100 s=subset(d, n>100) #### Descriptive statistics#### sum(d$n) xtabs(~d$Flipped) xtabs(~d$StudyID) range(xtabs(~d$StudyID)) mean(xtabs(~d$StudyID)) xtabs(~d$ReportedSignificant) xtabs(~d$Adjusted) xtabs(~d$ValencedEffect<= 0) xtabs(~d$ValencedEffect> 0) #Used to create labels on Fig1a fig1alabel1 = sum(xtabs(~d$ValencedEffect<= -0.2)["TRUE"]) fig1alabel2 = sum(xtabs(~d$ValencedEffect<= 0 & d$ValencedEffect> -0.2)["TRUE"]) fig1alabel3 = sum(xtabs(~d$ValencedEffect> 0 & d$ValencedEffect< 0.2)["TRUE"]) fig1alabel4 = sum(xtabs(~d$ValencedEffect>= 0.2)["TRUE"]) ###Histogram of correlations (fig 1a)#### fig1a=ggplot(d, aes(x=ValencedEffect)) + geom_histogram(binwidth=0.05, fill="grey70") + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black")) + coord_cartesian(xlim=c(-1, 1)) + xlab("Correlation with telomere measure") + ylab("Frequency") + geom_vline(xintercept=0) + geom_vline(xintercept=-0.2, linetype="dotted") + geom_vline(xintercept=0.2, linetype="dotted") + annotate("text", x=c(-0.5, -0.1, 0.1, 0.5), y=c(100, 100, 100, 100), label=c(fig1alabel1, fig1alabel2, fig1alabel3, fig1alabel4), size=3.5) fig1a win.metafile("fig1a.wmf", width=4, height=3, pointsize=14) fig1a dev.off() ### Whole dataset basic meta-analytic model #### m1=rma.mv(data=d, yi=ValencedEffect, V=CommonEffectVariance, random=~as.factor(number)|StudyID) summary(m1) ###Publication bias#### # Funnel plot (figure 1b) means=c( mean(d$ValencedEffect[d$ncat=="<=100"]), mean(d$ValencedEffect[d$ncat=="101-250"]), mean(d$ValencedEffect[d$ncat=="251-1000"]), mean(d$ValencedEffect[d$ncat==">1000"])) ns=c(50, 175, 675, 5000) fig1b=ggplot(d, aes(x=ValencedEffect, y=n)) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black")) + geom_point(shape=21, colour="black", fill="grey") + geom_vline(xintercept = 0) + scale_y_continuous(trans="sqrt", breaks=c(100, 500, 1000, 2500, 5000, 10000, 12500)) + ylab("Sample size") + xlab("Correlation with telomere measure") + annotate("point", x=means, y=ns, colour="red", size=4) + annotate("segment", x=means[1], xend=means[2], y=ns[1], yend=ns[2], colour="red", size=1.5) + annotate("segment", x=means[2], xend=means[3], y=ns[2], yend=ns[3], colour="red", size=1.5) + annotate("segment", x=means[3], xend=means[4], y=ns[3], yend=ns[4], colour="red", size=1.5) fig1b win.metafile("fig1b.wmf", width=4, height=3, pointsize=14) fig1b dev.off() #Quartiles of sample size quantile(d$n) xtabs(~d$ncat) # Separate meta-analyses by bin of sample size (figure 1c) s1=rma.mv(data=d, yi=ValencedEffect, V=CommonEffectVariance, random=~StudyID/number, subset=(ncat=="<=100")) s2=rma.mv(data=d, yi=ValencedEffect, V=CommonEffectVariance, random=~StudyID/number, subset=(ncat=="101-250")) s3=rma.mv(data=d, yi=ValencedEffect, V=CommonEffectVariance, random=~StudyID/number, subset=(ncat=="251-1000")) s4=rma.mv(data=d, yi=ValencedEffect, V=CommonEffectVariance, random=~StudyID/number, subset=(ncat==">1000")) d$ncat2=factor(d$ncat, levels=c("<=100", "101-250", "251-1000", ">1000")) ks=c(543, as.numeric(xtabs(~d$ncat2))) ms=c(length(unique(d$StudyID)), length(unique(d$StudyID[d$ncat=="<=100"])), length(unique(d$StudyID[d$ncat=="101-250"])), length(unique(d$StudyID[d$ncat=="251-1000"])), length(unique(d$StudyID[d$ncat==">1000"])) ) win.metafile("fig1c.wmf", width=8, height=4) forest.default(x=c(coef(m1), coef(s1), coef(s2), coef(s3), coef(s4)), vi=c(vcov(m1), vcov(s1), vcov(s2), vcov(s3), vcov(s4)), psize=0.5, slab=c("Whole dataset", "<=100", "101-250", "251-1000", ">1000"), ilab=cbind(ks, ms), ilab.xpos=c(0.05, 0.1), rows=c(6, 4:1), ylim=c(0,9), xlab="Correlation with telomere measure", alim=c(-0.3, 0.1)) text(-0.525, 7.5, "Sample size") text(0.05, 7.5, "k") text(0.1, 7.5, "m") text(0.25, 7.5, "r [95% CI]") dev.off() # Test whether correlations are different in different bins of sample size d$ncat=relevel(d$ncat, ref=">1000") m2=rma.mv(data=d, yi=ValencedEffect, V=CommonEffectVariance, random=~as.factor(number)|StudyID, mod=~ncat) summary(m2) ### Basic meta-analytic model for the reduced dataset#### length(unique(s$StudyID)) s1=rma.mv(data=s, yi=ValencedEffect, V=CommonEffectVariance, random=~as.factor(number)|StudyID) summary(s1) ###Analyses by category of exposure#### mcat = rma.mv(data=d, yi=ValencedEffect, V=CommonEffectVariance, random=~as.factor(number)|StudyID, mod=~AdversityOrEnrichmentType) summary(mcat) scat = rma.mv(data=s, yi=ValencedEffect, V=CommonEffectVariance, random=~as.factor(number)|StudyID, mod=~AdversityOrEnrichmentType) summary(scat) ### Figure 2 forest plot by broad category #### # Models on the reduced dataset by category disease.s = rma.mv(data=s, yi = ValencedEffect, V = CommonEffectVariance, random=~StudyID/number, subset = AdversityOrEnrichmentType=="disease") environmentalhazard.s = rma.mv(data=s, yi = ValencedEffect, V = CommonEffectVariance, random=~StudyID/number, subset = AdversityOrEnrichmentType=="environmentalhazard") nutrition.s = rma.mv(data=s, yi = ValencedEffect, V = CommonEffectVariance, random=~StudyID/number, subset = AdversityOrEnrichmentType=="nutrition") psychiatric.s = rma.mv(data=s, yi = ValencedEffect, V = CommonEffectVariance, random=~StudyID/number, subset = AdversityOrEnrichmentType=="psychiatric") smoking.s = rma.mv(data=s, yi = ValencedEffect, V = CommonEffectVariance, random=~StudyID/number, subset = AdversityOrEnrichmentType=="smoking") alcohol.s = rma.mv(data=s, yi = ValencedEffect, V = CommonEffectVariance, random=~StudyID/number, subset = AdversityOrEnrichmentType=="alcohol") sleep.s = rma.mv(data=s, yi = ValencedEffect, V = CommonEffectVariance, random=~StudyID/number, subset = AdversityOrEnrichmentType=="sleep") psychosocial.s = rma.mv(data=s, yi = ValencedEffect, V = CommonEffectVariance, random=~StudyID/number, subset = AdversityOrEnrichmentType=="psychosocial") physicalactivity.s=rma.mv(data=s, yi = ValencedEffect, V = CommonEffectVariance, random=~StudyID/number, subset = AdversityOrEnrichmentType=="physicalactivity") parentalcare.s=rma.mv(data=s, yi = ValencedEffect, V = CommonEffectVariance, random=~StudyID/number, subset = AdversityOrEnrichmentType=="parentalcare") socioeconomic.s=rma.mv(data=s, yi = ValencedEffect, V = CommonEffectVariance, random=~StudyID/number, subset = AdversityOrEnrichmentType=="socioeconomic") other.s=rma.mv(data=s, yi = ValencedEffect, V = CommonEffectVariance, random=~StudyID/number, subset = AdversityOrEnrichmentType=="other") # Models on the whole dataset by category disease.d = rma.mv(data=d, yi = ValencedEffect, V = CommonEffectVariance, random=~StudyID/number, subset = AdversityOrEnrichmentType=="disease") environmentalhazard.d = rma.mv(data=d, yi = ValencedEffect, V = CommonEffectVariance, random=~StudyID/number, subset = AdversityOrEnrichmentType=="environmentalhazard") nutrition.d = rma.mv(data=d, yi = ValencedEffect, V = CommonEffectVariance, random=~StudyID/number, subset = AdversityOrEnrichmentType=="nutrition") psychiatric.d = rma.mv(data=d, yi = ValencedEffect, V = CommonEffectVariance, random=~StudyID/number, subset = AdversityOrEnrichmentType=="psychiatric") smoking.d = rma.mv(data=d, yi = ValencedEffect, V = CommonEffectVariance, random=~StudyID/number, subset = AdversityOrEnrichmentType=="smoking") alcohol.d = rma.mv(data=d, yi = ValencedEffect, V = CommonEffectVariance, random=~StudyID/number, subset = AdversityOrEnrichmentType=="alcohol") sleep.d = rma.mv(data=d, yi = ValencedEffect, V = CommonEffectVariance, random=~StudyID/number, subset = AdversityOrEnrichmentType=="sleep") psychosocial.d = rma.mv(data=d, yi = ValencedEffect, V = CommonEffectVariance, random=~StudyID/number, subset = AdversityOrEnrichmentType=="psychosocial") physicalactivity.d=rma.mv(data=d, yi = ValencedEffect, V = CommonEffectVariance, random=~StudyID/number, subset = AdversityOrEnrichmentType=="physicalactivity") parentalcare.d=rma.mv(data=d, yi = ValencedEffect, V = CommonEffectVariance, random=~StudyID/number, subset = AdversityOrEnrichmentType=="parentalcare") socioeconomic.d=rma.mv(data=d, yi = ValencedEffect, V = CommonEffectVariance, random=~StudyID/number, subset = AdversityOrEnrichmentType=="socioeconomic") other.d=rma.mv(data=d, yi = ValencedEffect, V = CommonEffectVariance, random=~StudyID/number, subset = AdversityOrEnrichmentType=="other") other.d #Calculate median sample size for studies from each category in the whole dataset diseasen = median(subset(d, AdversityOrEnrichmentType=="disease")$n) environmentalhazardn = median(subset(d, AdversityOrEnrichmentType=="environmentalhazard")$n) nutritionn = median(subset(d, AdversityOrEnrichmentType=="nutrition")$n) psychiatricn = median(subset(d, AdversityOrEnrichmentType=="psychiatric")$n) smokingn = median(subset(d, AdversityOrEnrichmentType=="smoking")$n) alcoholn = median(subset(d, AdversityOrEnrichmentType=="alcohol")$n) sleepn = median(subset(d, AdversityOrEnrichmentType=="sleep")$n) psychosocialn = median(subset(d, AdversityOrEnrichmentType=="psychosocial")$n) physicalactivityn = median(subset(d, AdversityOrEnrichmentType=="physicalactivity")$n) parentalcaren = median(subset(d, AdversityOrEnrichmentType=="parentalcare")$n) socioeconomicn = median(subset(d, AdversityOrEnrichmentType=="socioeconomic")$n) othern = median(subset(d, AdversityOrEnrichmentType=="other")$n) #Calculate median sample size for studies from each category in the reduced dataset diseasens = median(subset(s, AdversityOrEnrichmentType=="disease")$n) environmentalhazardns = median(subset(s, AdversityOrEnrichmentType=="environmentalhazard")$n) nutritionns = median(subset(s, AdversityOrEnrichmentType=="nutrition")$n) psychiatricns = median(subset(s, AdversityOrEnrichmentType=="psychiatric")$n) smokingns = median(subset(s, AdversityOrEnrichmentType=="smoking")$n) alcoholns = median(subset(s, AdversityOrEnrichmentType=="alcohol")$n) sleepns = median(subset(s, AdversityOrEnrichmentType=="sleep")$n) psychosocialns = median(subset(s, AdversityOrEnrichmentType=="psychosocial")$n) physicalactivityns = median(subset(s, AdversityOrEnrichmentType=="physicalactivity")$n) parentalcarens = median(subset(s, AdversityOrEnrichmentType=="parentalcare")$n) socioeconomicns = median(subset(s, AdversityOrEnrichmentType=="socioeconomic")$n) otherns = median(subset(s, AdversityOrEnrichmentType=="other")$n) # Fig 2 forest plot names=c("Physical disease", " ", "Environmental hazard", " ", "Nutrition (poorer)", " ","Psychiatric illness", " ", "Smoking"," ", "Alcohol"," ", "Sleep (worse)", " ","Physical activity (less)"," ", "Psychosocial"," ", "Parental care (poorer)"," ", "Socioeconomic status (lower)"," ", "Other"," ") adversitytypeestimates= c(coef(disease.d), coef(disease.s), coef(environmentalhazard.d), coef(environmentalhazard.s), coef(nutrition.d), coef(nutrition.s), coef(psychiatric.d), coef(psychiatric.s), coef(smoking.d), coef(smoking.s), coef(alcohol.d), coef(alcohol.s), coef(sleep.d), coef(sleep.s), coef(physicalactivity.d), coef(physicalactivity.s), coef(psychosocial.d), coef(psychosocial.s), coef(parentalcare.d), coef(parentalcare.s), coef(socioeconomic.d), coef(socioeconomic.s), coef(other.d), coef(other.s)) adversitytypevariances= c(vcov(disease.d), vcov(disease.s), vcov(environmentalhazard.d), vcov(environmentalhazard.s), vcov(nutrition.d), vcov(nutrition.s), vcov(psychiatric.d), vcov(psychiatric.s), vcov(smoking.d), vcov(smoking.s), vcov(alcohol.d), vcov(alcohol.s), vcov(sleep.d), vcov(sleep.s), vcov(physicalactivity.d), vcov(physicalactivity.s), vcov(psychosocial.d), vcov(psychosocial.s), vcov(parentalcare.d), vcov(parentalcare.s), vcov(socioeconomic.d), vcov(socioeconomic.s), vcov(other.d), vcov(other.s)) studies=c(length(unique(d$StudyID[d$AdversityOrEnrichmentType=="disease"])), length(unique(s$StudyID[s$AdversityOrEnrichmentType=="disease"])), length(unique(d$StudyID[d$AdversityOrEnrichmentType=="environmentalhazard"])), length(unique(s$StudyID[s$AdversityOrEnrichmentType=="environmentalhazard"])), length(unique(d$StudyID[d$AdversityOrEnrichmentType=="nutrition"])), length(unique(s$StudyID[s$AdversityOrEnrichmentType=="nutrition"])), length(unique(d$StudyID[d$AdversityOrEnrichmentType=="psychiatric"])), length(unique(s$StudyID[s$AdversityOrEnrichmentType=="psychiatric"])), length(unique(d$StudyID[d$AdversityOrEnrichmentType=="smoking"])), length(unique(s$StudyID[s$AdversityOrEnrichmentType=="smoking"])), length(unique(d$StudyID[d$AdversityOrEnrichmentType=="alcohol"])), length(unique(s$StudyID[s$AdversityOrEnrichmentType=="alcohol"])), length(unique(d$StudyID[d$AdversityOrEnrichmentType=="sleep"])), length(unique(s$StudyID[s$AdversityOrEnrichmentType=="sleep"])), length(unique(d$StudyID[d$AdversityOrEnrichmentType=="physicalactivity"])), length(unique(s$StudyID[s$AdversityOrEnrichmentType=="physicalactivity"])), length(unique(d$StudyID[d$AdversityOrEnrichmentType=="psychosocial"])), length(unique(s$StudyID[s$AdversityOrEnrichmentType=="psychosocial"])), length(unique(d$StudyID[d$AdversityOrEnrichmentType=="parentalcare"])), length(unique(s$StudyID[s$AdversityOrEnrichmentType=="parentalcare"])), length(unique(d$StudyID[d$AdversityOrEnrichmentType=="socioeconomic"])), length(unique(s$StudyID[s$AdversityOrEnrichmentType=="socioeconomic"])), length(unique(d$StudyID[d$AdversityOrEnrichmentType=="other"])), length(unique(s$StudyID[s$AdversityOrEnrichmentType=="other"]))) ks=c(disease.d$k,disease.s$k, environmentalhazard.d$k,environmentalhazard.s$k, nutrition.d$k, nutrition.s$k, psychiatric.d$k,psychiatric.s$k, smoking.d$k, smoking.s$k, alcohol.d$k, alcohol.s$k, sleep.d$k, sleep.s$k, physicalactivity.d$k, physicalactivity.s$k, psychosocial.d$k, psychosocial.s$k, parentalcare.d$k, parentalcare.s$k,socioeconomic.d$k,socioeconomic.s$k, other.d$k, other.s$k) ns=c(diseasen,diseasens, environmentalhazardn,environmentalhazardns, nutritionn, nutritionns, psychiatricn, psychiatricns, smokingn, smokingns, alcoholn, alcoholns, sleepn, sleepns, physicalactivityn, physicalactivityns, psychosocialn, psychosocialns, parentalcaren, parentalcarens, socioeconomicn,socioeconomicns, othern, otherns) fig2dat=data.frame(adversitytypeestimates, adversitytypevariances, names, ks, studies, ns) png("fig2.png", res=300, units="in", width=10, height=8) forest(x=fig2dat$adversitytypeestimates, vi=fig2dat$adversitytypevariances, slab=fig2dat$names, xlab="Correlation with telomere measure", psize=1, col=rep(c("black", "grey50"), times=12), ilab=cbind(fig2dat$ns, fig2dat$ks, fig2dat$studies), ilab.xpos=c(-0.6, 0.3, 0.4)) text(-1.125, 25.5, "Type of exposure") text(-0.6, 25.5, "Median sample size") text(0.3, 25.5, "k") text(0.4, 25.5, "m") text(0.7, 25.5, "r [95% CI]") for (i in (seq(2.5, 22.5, by=2))) {abline(h=i, lty=3, col="grey")} dev.off() ###Analysis by fine categories of exposure #### # Whole dataset d$FineAdversityType2=relevel(d$FineAdversityType, ref="Smoking") m.fine=rma.mv(data=d, yi=ValencedEffect, V=CommonEffectVariance, random=~as.factor(number)|StudyID, mods=~FineAdversityType2) m.fine # Reduced dataset s$FineAdversityType2=relevel(s$FineAdversityType, ref="Smoking") m.fine.r=rma.mv(data=s, yi=ValencedEffect, V=CommonEffectVariance, random=~as.factor(number)|StudyID, mods=~FineAdversityType2) m.fine.r ### Figure 3: The mother of all forest plots#### new.order=c(3, 5, 10, 16, 23, 29, 12, 4, 8, 13, 14, 15, 18, 25, 33, 35, 2, 9, 26, 27, 30, 1, 28, 24, 6, 19, 31, 32, 34, 36, 22, 7, 11, 17, 20, 21) d$FineAdversityType2=factor(d$FineAdversityType, levels(d$FineAdversityType)[new.order]) s=subset(d, n>100) i=0 names=NULL coeffs=NULL variances=NULL ks=NULL ms=NULL adj=NULL ns=0 for (current in levels(d$FineAdversityType2)) { i = i+1 names[i]=current m=rma.mv(data=d, yi=ValencedEffect, V=CommonEffectVariance, random=~StudyID/number, subset=(FineAdversityType2==current)) coeffs[i]=coef(m) variances[i]=vcov(m) ks[i]=m$k ms[i]=length(unique(d$StudyID[d$FineAdversityType2==current])) ns[i]=median(subset(d, FineAdversityType2==current)$n) i = i+1 names[i]=" " if (current=="Mood (lower)"|current=="Vitamins (lower)"|current=="Caregiving") {coeffs[i]=0 variances[i]=0 ks[i]=0 ms[i]=0 ns[i]=0 } if (current!="Mood (lower)" & current!="Vitamins (lower)" & current!="Caregiving") { m=rma.mv(data=d, yi=ValencedEffect, V=CommonEffectVariance, random=~StudyID/number, subset=(FineAdversityType2==current & ncat!="<=100")) coeffs[i]=coef(m) variances[i]=vcov(m) ks[i]=m$k ms[i]=length(unique(d$StudyID[d$FineAdversityType2==current & d$ncat!="<=100"])) ns[i]=median(subset(s, FineAdversityType2==current)$n) } } psizes=rep(1, times=72) psizes[which(coeffs==0)]=0 png("figure3.png", res=300, width=10, height=13, units="in") forest.default(x = coeffs, vi=variances, slab=names, psize=psizes, col=rep(c("black", "grey50"), times=length(names)/2), alim=c(-0.6, 0.2), xlim=c(-1, 1), xlab="Correlation with telomere measure", ilab=cbind(ns, ks, ms), ilab.xpos=c(0.3, 0.5, 0.6), cex=1 ) for (i in (seq(2.5, 71.5, by=2))) {abline(h=i, lty=3, col="grey")} text(-0.8, 74, "Exposure fine category") text(0.875, 74, "r [95% CI]") text(0.5, 74, "k") text(0.6, 74, "m") text(0.3, 74, "Median sample size") dev.off() win.metafile(filename = "figure3.wmf", width=10, height=15) forest.default(x = coeffs, vi=variances, slab=names, psize=psizes, col=rep(c("black", "grey50"), times=length(names)/2), alim=c(-0.6, 0.2), xlim=c(-1, 1), xlab="Correlation with telomere measure", ilab=cbind(ns, ks, ms), ilab.xpos=c(0.3, 0.5, 0.6), cex=1 ) for (i in (seq(2.5, 71.5, by=2))) {abline(h=i, lty=3, col="grey")} text(-0.8, 74, "Exposure fine category") text(0.875, 74, "r [95% CI]") text(0.5, 74, "k") text(0.6, 74, "m") text(0.3, 74, "Median sample size") dev.off() ###Other moderators (table 2)##### # m models use full dataset, s models use reduced dataset # Longitudinal xtabs(~d$Longitudinal) mo3=rma.mv(data=d, yi=ValencedEffect, V=CommonEffectVariance, random=~StudyID/number, mod=~Longitudinal) summary(mo3) so3=rma.mv(data=s, yi=ValencedEffect, V=CommonEffectVariance, random=~StudyID/number, mod=~Longitudinal) summary(so3) # Experimental xtabs(~d$Experimental) mo4=rma.mv(data=d, yi=ValencedEffect, V=CommonEffectVariance, random=~StudyID/number, mod=~Experimental) summary(mo4) so4=rma.mv(data=s, yi=ValencedEffect, V=CommonEffectVariance, random=~StudyID/number, mod=~Experimental) summary(so4) # Lifestage exposure xtabs(~d$Lifestage) mo5=rma.mv(data=d, yi=ValencedEffect, V=CommonEffectVariance, random=~StudyID/number, mod=~Lifestage) summary(mo5) so5=rma.mv(data=s, yi=ValencedEffect, V=CommonEffectVariance, random=~StudyID/number, mod=~Lifestage) summary(so5) # Now try without the lifestage unreporteds mo5a=rma.mv(data=d, yi=ValencedEffect, V=CommonEffectVariance, random=~StudyID/number, mod=~Lifestage, subset=(Lifestage!="unreported")) summary(mo5a) so5a=rma.mv(data=s, yi=ValencedEffect, V=CommonEffectVariance, random=~StudyID/number, mod=~Lifestage, subset=(Lifestage!="unreported")) summary(so5a) # Lifestage measurement xtabs(~d$LifestageAtLastTL) mo6=rma.mv(data=d, yi=ValencedEffect, V=CommonEffectVariance, random=~StudyID/number, mod=~LifestageAtLastTL) summary(mo6) so6=rma.mv(data=s, yi=ValencedEffect, V=CommonEffectVariance, random=~StudyID/number, mod=~LifestageAtLastTL) summary(so6) # Now try without the unreporteds mo6a=rma.mv(data=d, yi=ValencedEffect, V=CommonEffectVariance, random=~StudyID/number, mod=~LifestageAtLastTL, subset=(LifestageAtLastTL!="unreported")) summary(mo6a) so6a=rma.mv(data=s, yi=ValencedEffect, V=CommonEffectVariance, random=~StudyID/number, mod=~LifestageAtLastTL, subset=(LifestageAtLastTL!="unreported")) summary(so6a) xtabs(~s$LifestageAtLastTL) # There's only one study in the reduced dataset with lifestage at last TL = embryonic # Tissue type xtabs(~d$TissueType) d$TissueType=relevel(d$TissueType, ref="wbcs") mo7=rma.mv(data=d, yi=ValencedEffect, V=CommonEffectVariance, random=~StudyID/number, mod=~TissueType) summary(mo7) s$TissueType=relevel(s$TissueType, ref="wbcs") so7=rma.mv(data=s, yi=ValencedEffect, V=CommonEffectVariance, random=~StudyID/number, mod=~TissueType) summary(so7) xtabs(~s$TissueType) # Technique describeBy(d$n, d$TechniqueType) describeBy(s$n, s$TechniqueType) d$TechniqueType=relevel(d$TechniqueType, ref="qpcr") s$TechniqueType=relevel(s$TechniqueType, ref="qpcr") # Full dataset mo8=rma.mv(data=d, yi=ValencedEffect, V=CommonEffectVariance, random=~StudyID/number, mod=~TechniqueType) summary(mo8) # Reduced dataset so8=rma.mv(data=s, yi=ValencedEffect, V=CommonEffectVariance, random=~StudyID/number, mod=~TechniqueType) summary(so8) # Full dataset controlling for sample size x08=rma.mv(data=d, yi=ValencedEffect, V=CommonEffectVariance, random=~StudyID/number, mod=~TechniqueType + sqrt(n), btt=2:4) summary(x08) # The moderating effect of technique disappears once you control for sample size # Separate metaanalyses by technique qpcr=rma.mv(data=d, yi=ValencedEffect, V=CommonEffectVariance, random=~StudyID/number, subset=(TechniqueType=="qpcr")) summary(qpcr) southern=rma.mv(data=d, yi=ValencedEffect, V=CommonEffectVariance, random=~StudyID/number, subset=(TechniqueType=="southernblot")) summary(southern) fish=rma.mv(data=d, yi=ValencedEffect, V=CommonEffectVariance, random=~StudyID/number, subset=(TechniqueType=="fish")) summary(fish) telseq=rma.mv(data=d, yi=ValencedEffect, V=CommonEffectVariance, random=~number, subset=(TechniqueType=="telseq")) summary(telseq) # Extra figure showing the central estimates for each technique type (not used in paper) fig4dat=data.frame( c("qPCR", "Southern blot", "FISH", "Telseq"), c(coef(qpcr), coef(southern), coef(fish), coef(telseq)), c(vcov(qpcr), vcov(southern), vcov(fish), vcov(telseq)), c(qpcr$k, southern$k, fish$k, telseq$k), c( length(unique(d$StudyID[d$TechniqueType=="qpcr"])), length(unique(d$StudyID[d$TechniqueType=="southernblot"])), length(unique(d$StudyID[d$TechniqueType=="fish"])), length(unique(d$StudyID[d$TechniqueType=="telseq"])) )) names(fig4dat)=c("Technique", "Estimate", "V", "k", "m") forest.default(x=fig4dat$Estimate, vi=fig4dat$V, slab=fig4dat$Technique, xlab="Correlation with telomere measure", ilab=paste(fig4dat$k,"(",fig4dat$m, ")"), ilab.xpos = 0.175, psize=0.5, alim=c(-0.6, 0.2), xlim=c(-.8, 0.6)) text(-0.7, 5.5, "Technique") text(0.175, 5.5, "k (m)") text(0.45, 5.5, "r [95% CI]") # Sex mo9=rma.mv(data=d, yi=ValencedEffect, V=CommonEffectVariance, mod=~Sex, random=~1|StudyID/number) summary(mo9) so9=rma.mv(data=s, yi=ValencedEffect, V=CommonEffectVariance, mod=~Sex, random=~1|StudyID/number) summary(so9) ###Discussion#### length(d$n[d$n>359])/length(d$n) length(d$n[d$n>1007])/length(d$n) summarise(group_by(d, Longitudinal), followup=mean(as.numeric(as.character(AgeAtLastTL))-as.numeric(as.character(AgeAtFirstTL)), na.rm=T), sd=sd(as.numeric(as.character(AgeAtLastTL))-as.numeric(as.character(AgeAtFirstTL)), na.rm=T)) ####Supplementary analyses (as reported in supplementary document) ####### ###Supplement section 1#### # Compare estimates incorporating and ignoring the multilevel structureof the data m1=rma.mv(data=d, yi=ValencedEffect, V=CommonEffectVariance, random=~factor(number)|StudyID) summary(m1) x1=rma(data=d, yi=ValencedEffect, vi=CommonEffectVariance) summary(x1) # Studies with more associations also have weaker ones require(dplyr) meansum=summarise(group_by(d, StudyID), nassociations=length(ValencedEffect), meanassociation=mean(ValencedEffect)) cor.test(meansum$nassociations, abs(meansum$meanassociation)) # Function to create a flat (not multilevel) version of the dataset by randomly choosing one association per study flat=function(d) {i = levels(d$StudyID)[1] current.subset=subset(d, StudyID==i) nrows=dim(current.subset)[1] f=current.subset[sample(nrows, 1), ] for (i in levels(d$StudyID)[2:length(levels(d$StudyID))]) { current.subset=subset(d, StudyID==i) nrows=dim(current.subset)[1] chosen.row=current.subset[sample(nrows, 1), ] f=rbind(f, chosen.row) } return(f) } # Comparison of results from flat datasets, multilevel analysis, and ignoring multilevel structure # Note you will get slightly different results each time due to sampling central=NULL ses=NULL central[1]=m1$b ses[1]=m1$se for (i in 2:11) { fi=flat(d) xi=rma(fi$ValencedEffect, vi=fi$CommonEffectVariance) central[i]=xi$b ses[i]=xi$se } central[12]=rma(data=d, yi=ValencedEffect, vi=CommonEffectVariance)$b ses[12]=rma(data=d, yi=ValencedEffect, vi=CommonEffectVariance)$se # Figure S1 png("figures1.png", res=300, units="in", width=6, height = 4) forest.default(x=central, sei=ses, alim=c(-0.25, 0), xlim=c(-0.5, 0.25), psize=0.75, slab=c("Multilevel model", "Flat 1", "Flat 2", "Flat 3", "Flat 4", "Flat 5", "Flat 6", "Flat 7", "Flat 8", "Flat 9", "Flat 10", "Ignoring multilevel"), xlab="Correlation with telomere measure") dev.off() ###Supplement section 2: Mixture model to identify outliers#### # NB This code takes several hours to run require(metaplus) # This next code samples the flat dataset ten times and fits both a normal and a mixture model to each normal=NULL mixture=NULL BIC.normal=NULL BIC.mixture=NULL tau=NULL tau.out=NULL out.p=rep(0, times=length(unique(d$StudyID))) for (i in 1:10) { fi=flat(d) mp1 = metaplus(data=fi, yi=ValencedEffect, sei=sqrt(CommonEffectVariance)) mp2 = metaplus(data=fi, yi=ValencedEffect, sei=sqrt(CommonEffectVariance), random="mixture") normal[i]=mp1$results[1, 1] mixture[i]=mp2$results[1, 1] tau[i]=mp2$results[2, 1] tau.out[i]=mp2$results[3, 1] BIC.normal[i]=summary(mp1)$fitstats$BIC BIC.mixture[i]=summary(mp2)$fitstats$BIC out.p=out.p+outlierProbs(mp2)$'outlier.prob' print(i) } out.p=out.p/10 # Write the main results of ecah run to a file metaplus.output=data.frame(normal, mixture, BIC.normal, BIC.mixture, tau, tau.out) write.csv(metaplus.output, "metaplus.output.csv") # Write the outlier probability for each study to a file outliers = data.frame(fi$StudyID, out.p) write.csv(outliers, "outlier.probabilities.csv") # Fit for standard and mixture models mean(BIC.mixture) mean(BIC.normal) # Dispersion of outlier and non-outlier studies mean(sqrt(tau)) mean(sqrt(tau.out)) # Estimates from mixture and standard model mean(mixture) mean(normal) # This tells you which studies are consistently outliers fi$StudyID[out.p>0.9] # And their sample sizes fi$n[out.p>0.9] median(fi$n[out.p>0.9]) # And their association strengths fi$ValencedEffect[out.p>0.9] ###Supplement section 3: Publication bias model using weightr#### # NB results will not be numerically identical every time due to sampling require(weightr) # This section samples ten times from the whole dataset and runs the weightr publication bias analysis on each unadjustedcentral=NULL central=NULL antitrend=NULL protrend=NULL sig=NULL for (i in 1:10) { fi=flat(d) wfi=weightfunct(effect = fi$ValencedEffect, v = fi$CommonEffectVariance, table=T, steps=c(0.025, 0.5, 0.975)) unadjustedcentral[i]=wfi[[1]]$par[2] central[i]=wfi[[2]]$par[2] antitrend[i]=wfi[[2]]$par[3] protrend[i]=wfi[[2]]$par[4] sig[i]=wfi[[2]]$par[5] } # The next two lines give you the mean association in unadjusted vs. adjusted model mean(unadjustedcentral) mean(central) # The next lines give you the probability of observation of different kinds of association # Reference category is a significant association in the opposite direction to general expectation mean(antitrend) mean(protrend) mean(sig) # Figure S2: Showing the funnel plot again but labelled by whether the results were reported significant by the original authors means=c( mean(d$ValencedEffect[d$ncat=="<=100"]), mean(d$ValencedEffect[d$ncat=="101-250"]), mean(d$ValencedEffect[d$ncat=="251-1000"]), mean(d$ValencedEffect[d$ncat==">1000"])) ns=c(50, 175, 675, 5000) figS2=ggplot(d, aes(x=ValencedEffect, y=n)) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), legend.position=c(0.8, 0.8)) + annotate("point", x=means, y=ns, colour="grey", size=4) + annotate("segment", x=means[1], xend=means[2], y=ns[1], yend=ns[2], colour="grey", size=1.5) + annotate("segment", x=means[2], xend=means[3], y=ns[2], yend=ns[3], colour="grey", size=1.5) + annotate("segment", x=means[3], xend=means[4], y=ns[3], yend=ns[4], colour="grey", size=1.5) + geom_point(shape=21, colour="black", aes(fill=ReportedSignificant)) + geom_vline(xintercept = 0) + scale_y_continuous(trans="sqrt", breaks=c(100, 500, 1000, 2500, 5000, 10000, 12500)) + ylab("Sample size") + xlab("Correlation with telomere measure") + scale_fill_manual(name="Reported significant", labels=c("No", "Yes", "Unreported"), values=c("red", "green", "white")) figS2 png("figureS2.png", res=300, width=6, height=4, units="in") figS2 dev.off()