Commit 8e54ba17 authored by Corson N. Areshenkoff's avatar Corson N. Areshenkoff

Initial commit

parent 9b49e157
# 2017_smart_metaAnalysis # smartMetaAnalysis
Data and code for replicating the meta-analysis reported in Smart et al. (2017).
The file 'data/masterSpreadsheet.csv' reports raw outcomes reported by the studies
included in the meta-analysis. The meta-analysis models are located in the 'models'
folder.
Code for replicating the analyses reported in Smart, et al. 2017. The script 'runMetaAnalysis.R' runs the complete analysis, from the importing and preparation
\ No newline at end of file of the raw data in 'data/masterSpreadsheet.csv', to the plotting of the fitted model. More
details are found within the '.R' files.
Currently, the code only replicates the cognitive outcome model. Code for the full model
is forthcoming.
## References
Smart, C. M., Karr, J. E., Areshenkoff, C. N., Rabin, L. A., Hudon, C., Gates, N., ... & Hampel, H. (2017). Non-pharmacologic interventions for older adults with subjective cognitive decline: systematic review, meta-analysis, and preliminary recommendations. Neuropsychology Review, 1-13.
library(plyr)
# Load data spreadsheet
raw <- read.table('data/masterSpreadsheet.csv', sep = ',',
header = T, na.strings = 'NR')
# Study info frame
studyInfo <- ddply(raw, c("Author", "Condition"), summarise,
Year = Year[1],
N = Sample.Size[1],
Age = Age[1],
Group.or.Individual = Group.or.Individual[1],
Control.Type = Control[1],
Percent.Male = Percent.Male[1],
Percent.White = Percent.White[1],
MMSE = MMSE[1],
PEDro = PEDro[1],
Treatment.Length = Treatment.Length[1],
Treatment.Frequency = Treatment.Frequency[1],
Session.Duration = Tx.Sesson.Duration[1])
# Compute mean differences
# Note that some studies report both pre- and post-test means,
# while others report pre-test means and mean differences.
# For the former, we compute the differences here.
for (i in 1:dim(raw)[1]){
if (is.na(raw$Mean.Change[i])){
raw$Mean.Change[i] <- raw$Mean.Post[i] - raw$Mean.Pre[i]
}
}
# Comput effect sizes and variances
# For each study, we compute the Morris effect size using
# the control as a baseline.
raw$Effect.Size <- NA
raw$ES.Variance <- NA
numStudies <- length(levels(raw$Author))
studies <- levels(raw$Author)
for (i in 1:numStudies){
dummyStudy <- subset(raw, Author == studies[i])
# Split by condition
numConditions <- length(levels(factor(dummyStudy$Condition))) - 1
treatmentIdx <- which(levels(factor(dummyStudy$Condition)) != 'Control')
conditions <- levels(factor(dummyStudy$Condition))[treatmentIdx]
dummyControl <- subset(dummyStudy, Condition == 'Control')
for (j in 1:numConditions){
# Compute effect size
dummyCondition <- subset(dummyStudy, Condition == conditions[j])
nT <- dummyCondition$Sample.Size
nC <- dummyControl$Sample.Size
cP <- 1 - 3 / ( 4*(nT + nC - 2) - 1 )
SD <- sqrt( ( (nT-1)*dummyCondition$SD.Pre^2 + (nC-1)*dummyControl$SD.Pre^2 ) /
(nT+nC-2))
ES <- cP * ( ( dummyCondition$Mean.Change - dummyControl$Mean.Change ) / SD)
raw$Effect.Size[raw$Author == studies[i] & raw$Condition == conditions[j]] <- ES
# Compute variance
raw$ES.Variance[raw$Author == studies[i] & raw$Condition == conditions[j]] <-
(nT+nC)/(nT*nC) + (ES^2)/(2*(nT + nC - 2))
}
}
# Remove the control studies
maData <- subset(raw, Condition != 'Control')
# Flip effects sizes to that positive is better
maData$Effect.Size <- maData$Effect.Size * maData$Higher.Score.Direction
# Save
write.table(studyInfo, file = 'data/studyInfo.csv',
sep = ',', row.names = F)
write.table(maData, file = 'data/maData.csv',
sep = ',', row.names = F)
# Load sample
load('data/cognitiveModelSamples')
samples = rstan::extract(fit)
# Labels
plotFrame = data.frame(Study = as.character(unique(fitData$Author)),
Mean = colMeans(samples$theta),
Lower = apply(samples$theta, 2, function(x) quantile(x, .025)),
Upper = apply(samples$theta, 2, function(x) quantile(x, .975)))
# We plot the studies in order of their mean effect size,
# so we get the ordering here
idx = order(plotFrame$Mean)
plotFrame = plotFrame[order(plotFrame$Mean),]
## The rest of the columns in the table.
tabletext <- cbind(c("Study", "\n", as.character(plotFrame$Study), '\n', 'Overall\n', ''))
pdf("plots/cognitive_forestPlot.pdf",
width = 6, height = 5)
forestplot(labeltext = tabletext, graph.pos=2,
mean=c(NA, NA, c(plotFrame$Mean, NA, mean(samples$mu), NA)),
lower = c(NA, NA, plotFrame$Lower, NA, quantile(samples$mu, .025), NA),
upper=c(NA, NA, plotFrame$Upper, NA, quantile(samples$mu, .975), NA),
title = "",
xlab = "Effect Size",
hrzl_lines = list("2" = gpar(lwd=2, col="#99999922")),
col=fpColors(box = "black", lines = "black", zero = "gray50"),
zero = 0, cex = 0.9, lineheight = "auto",
colgap = unit(6,"mm"), cex.lab = 1.5,
lwd.ci = 2, ci.vertices = F, grid = T,
txt_gp = fpTxtGp(ticks = gpar(cex=1),
xlab = gpar(cex=1)))
dev.off()
maData = read.table('data/maData.csv', sep = ',', header = T)
# Load data
fitData <- maData[maData$Method == 'Cognitive' &
maData$Condition == 'Cognitive',]
fitData$Author <- factor(fitData$Author,
levels = unique(as.character(fitData$Author)))
# Prepare data for Stan
stanData <- list(N = dim(fitData)[1],
K = length(levels(fitData$Author)),
Y = fitData$Effect.Size,
V = fitData$ES.Variance,
Study = as.numeric(fitData$Author))
# Fit model
fit = stan(file = 'models/maModel.stan',
data = stanData)
save(fit, file = 'data/cognitiveModelSamples')
\ No newline at end of file
This diff is collapsed.
data {
int<lower=0> N; // Number of effect sizes
int<lower=0> K; // Number of Studies
vector[N] Y; // Effect sizes
vector<lower=0>[N] V; // Variance of ES
int<lower=0> Study[N]; // Study indicator
}
parameters {
vector[K] theta_dummy; // Study effects
vector[N] d_dummy; // True effects
vector<lower=0>[K] sigma; // Study level variance
real mu; // Mean effect size
real<lower=0> tau; // Between study variance
}
transformed parameters {
vector[K] theta;
vector[N] d;
for (i in 1:K){
theta[i] = mu + theta_dummy[i] * tau;
}
for (i in 1:N){
d[i] = theta[Study[i]] + d_dummy[i] * sigma[Study[i]];
}
}
model {
// Model
for (i in 1:N){
Y[i] ~ normal(d[i], sqrt(V[i]));
}
d_dummy ~ normal(0,1);
theta_dummy ~ normal(0,1);
// Prior
mu ~ normal(0, 1);
tau ~ student_t(3, 0, 0.5);
sigma ~ inv_gamma(3, 0.5);
}
\ No newline at end of file
Github doesn't like empty folders
library(rstan)
library(forestplot)
library(plyr)
# Load raw data and calculate effect sizes
# The script generates two spreadsheets in the 'data/' folder
# maData.csv: Contains effect sizes and associated information
# for each outcome in each study.
# studyInfo.csv: Contains information about each study
source('cleanData.R')
rm(list = ls())
# Fit the model for the cognitive outcomes
# This script imports maData.csv and fits the meta-analysis
# model described in Smart et al. (2017) using the model code
# in 'models/maModel.stan'. The fitted model is saved in
# 'data/cognitiveModelSamples'
source('cognitiveModelFit.R')
# Finally, create a forest plot of the posterior means and
# 95% posterior intervals. The plot is saved as
# 'plots/cognitive_forestPlot.pdf'
source('cognitiveForestPlot.R')
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment