Compare Methods for Marine Meta-Analysis

method_comparison_LRR1.R

#######################################################################################################################################
#
# Compare Different Types of Meta-Analysis Method
# for calculating an effect size from LRR1
#
# Jarrett Byrnes
# Last Updated Dec 9,2012 11am
#
# Changelog
#######################################################################################################################################

## Load data and Code
marine <- read.csv("./BEF_MetaMaster marine merged_20121110.csv", na.strings=c(".", " "))

#############################################
##DATA FILTERING
marine=marine[-c(11,17,69,123,155),]

#Remove rows where Ycat is not consumption/production/biogeo_flux
marine=marine[marine$Ycat=="Consumption" | marine$Ycat=="Production" | marine$Ycat=="Biogeo_flux",]

#Convert certain columns to numeric?
#marine[,c(28:38,73:ncol(marine))]=apply(marine[,c(28:38,73:ncol(marine))],2,function(x) as.numeric(as.character(x)))

marine <- marine[-155,] #Emmett's study - entry 160, sty.id.marine=29, expt.id.marine=84 with bizarre low LRR
#############################################

#SUBSET TO JUST CONSUMPTION
cons <- subset(marine, marine$Ycat=="Consumption")
cons <- cons[-which(is.na(cons$LRR1)),]



#LOAD LIBRARIES
#library(devtools)#JON RUN THIS
#install_github("robustmeta", username="jebyrnes") #JON RUN THIS
library(metafor)
library(robustmeta)
library(nlme)
library(ggplot2)

#RUN A NUMBER OF MODELS UNDER DIFFERENT ASSUMPTIONS
nogroup <- rma(LRR1, VLRR1, data=cons)
group <- rrma(LRR1 ~ 1, study_id = Study, data=cons, var_eff = VLRR1, rho=0.2)  
unweighted <- lme(LRR1 ~ 1, random =~1|Study, data=cons)
cons$invVLRR1 <- 1/(cons$VLRR1+0.001)
cons$invN <- 1/(cons$NSmax)
ran_varweighted <- lme(LRR1 ~ 1, random =~1|Study, data=cons, weights=~invVLRR1)
ran_sampweighted <- lme(LRR1 ~ 1, random =~1|Study, data=cons, weights=~invN)

#EXTRACT MEANS AND CIS
coef_df <- data.frame(Type = c("Var Weighted Metafor", "Group/Var Weighted", "Mixed Unweighted", "Mixed Var Weighted", "Mixed N Weighted"),
         Estimate = c(predict(nogroup)$pred,
               predict(group)[1],
               fixef(unweighted)[1],
               fixef(ran_varweighted)[1],
               fixef(ran_sampweighted)[1]),
           
           CI.LB = c(predict(nogroup)$ci.lb,
                     group$est$ci.lb,
                     intervals(unweighted)$fixed[1],
                     intervals(ran_varweighted)$fixed[1],
                     intervals(ran_sampweighted)$fixed[1]),
           CI.UB = c(predict(nogroup)$ci.ub,
                     group$est$ci.ub,
                     intervals(unweighted)$fixed[3],
                     intervals(ran_varweighted)$fixed[3],
                     intervals(ran_sampweighted)$fixed[3])
           

)


#PLOT THE COMPARISON
library(ggplot2)

ggplot(data=coef_df, mapping=aes(x=Type, y=Estimate, ymin=CI.LB, ymax=CI.UB, color=Type)) +
  geom_pointrange(size=1.5) + theme_bw() +
  geom_hline(h=0, lty=2, col="red") +
  coord_flip()