# Janaina Silva Fortirer, 2018

############################################################################################################################
#FUNCAO CALCULA MEDIA PADRONIZADA POR ESTUDO, INTERVALO DE CONFIANCA E CRIA UM GRAFICO FOREST PLOT PARA META ANALISE
metaR<- function(metadata){

############################################################################################################################  
#VERIFICA OS PARAMETROS

# verifica se a classe dos valores da primeira coluna eh 'factor' 
if (class(metadata[,1]) == "factor")
{  
  metadata[,1] <- as.character(metadata[,1])# se sim, modifica a classe para caracter. obs. eh necessario que o nome dos autores estejam em caracter para serem utilizados na legenda do grafico
}
# verifica se a quantidade de linhas do data.frame eh menor que 2
if ((sum(row(metadata))) < 2)
{
  stop("To calculate effect size by meta-analysis must have at least 2 studies, with least 2 'row'") # para e exibe mensagem de erro ao usuario
}
# verifica se existe algum NA nos dados
if (sum(rowSums(is.na(metadata))) > 0)
{
 stop("The 'metadata' must not 'NAs'") # se sim, exibe mensagem de erro ao usuario
}
# verifica se argumento metadata eh do tipo data.frame
if(!is.data.frame(metadata)) 
{
 stop ("The argument metadata must be a 'data.frame'") # Para funcao e exibe mensagem de erro.
}
# Verifica se o nome na linha 1 e coluna 1 foi inserido corretamente
if (colnames(metadata[1]) != "Author") 
{
 stop ("The name in row 1 and column 1 must be 'Author'") # Para funcao e exibe mensagem de erro.
}
# Verifica se o nome na linha 1 e coluna 2 foi inserido corretamente
if (colnames(metadata[2]) != "Year") 
{
 stop ("The name in row 1 and column 2 must be 'Year'") # Para funcao e exibe mensagem de erro.
}
# Verifica se o nome na linha 1 e coluna 3 foi inserido corretamente
if (colnames(metadata[3]) != "Ncontrol") 
{
  stop ("The name in row 1 and column 1 must be 'Ncontrol'") # Para funcao e exibe mensagem de erro.
}
# Verifica se o nome na linha 1 e coluna 4 foi inserido corretamente
if (colnames(metadata[4]) != "Ntreated") 
{
  stop ("The name in row 1 and column 1 must be 'Ntreated'") # Para funcao e exibe mensagem de erro.
}
# Verifica se o nome na linha 1 e coluna 5 foi inserido corretamente
if (colnames(metadata[5]) != "MeanControl") 
{
  stop ("The name in row 1 and column 5 must be 'MeanControl'") # Para funcao e exibe mensagem de erro.
}
# Verifica se o nome na linha 1 e coluna 6 foi inserido corretamente
if (colnames(metadata[6]) != "MeanTreated") 
{
  stop ("The name in row 1 and column 6 must be 'MeanTreated'") # Para funcao e exibe mensagem de erro.
}
# Verifica se o nome na linha 1 e coluna 7 foi inserido corretamente
if (colnames(metadata[7]) != "Sdcontrol")
{
  stop ("The name in row 1 and column 7 must be 'SDcontrol'") # Para funcao e exibe mensagem de erro.
}
# Verifica se o nome na linha 1 e coluna 7 foi inserido corretamente
if (colnames(metadata[8]) != "Sdtreated") 
{
  stop ("The name in row 1 and column 7 must be 'SDtreated'") # Para funcao e exibe mensagem de erro.
}

############################################################################################################################
#INICIO DO CALCULO DA ESTATISTICA DE INTERESSE (EFFECT SIZE - 'd')
  
# calculo df (grau de liberdade) para calcular effect size 'd'
metadata$df <- with(metadata, (Ntreated + Ncontrol -2))
  
# calculo do fator de correcao (J), que sao sempre menor que 1.0 
# J geralmente tem valores proximos de 1.0, a menos que o df (grau de liberdade) seja < 10 (Hedges, 1981).
# formula para fator de correcao "J" (Hedges and Olkin, 1985)
metadata$J <- with(metadata, (1 - 3/(4 * (df) - 1))) 

# calculo do efeito de interesse 'd', que eh a diferenca entre as medias padronizada do resultado para cada estudo,
# formula de Hedges and Olkin, 1985. 
metadata$d <- with(metadata, (MeanTreated - MeanControl)/sqrt(((Ntreated - 1) * Sdtreated^2 + (Ncontrol - 1) * Sdcontrol^2)/(Ntreated + Ncontrol - 2)) *J ) 
  
# calculo da variancia de 'd' (formula: Hedges and Olkin, 1985)
metadata$var.d <- with(metadata, ((Ntreated + Ncontrol/Ntreated * Ncontrol) + (d^2 /2 * (Ntreated + Ncontrol))))
  
# o calculo da variancia tem importancia para a estatistica do efeito de interesse, 
# porem como o 'd' teve uma metrica quadratica, a interpretacao nao eh muito intuitiva
# por isso eh recomendado por (Borenstein et al. 2009) usar o erro padrao que esta na mesma escala que o 'd'
  
# Entao, iremos calcular do erro padrao de 'd'
metadata$SE.d <- with(metadata, sqrt(var.d))
  
###########################################################################################################################
# CALCULO DO INTERVALO DE CONFIANCA

# Nestas equacoes 1.96 eh o valor 'Z' correspondente aos limites de confianca de 95%, 
# isso permite o erro de 2.5% em cada extremidade da distribuicao (Borenstein et al. 2009).
  
# calculo do intervalo de confianca lower (CI - 95% )
metadata$CI.lower <- with (metadata, (d - 1.96 * SE.d))
  
  
# calculando o intervalo de confianca upper (CI - 95% )
metadata$CI.upper <- with(metadata, (d + 1.96 * SE.d))
  
############################################################################################################################
# CALCULO DA ESTATISTICA 'Z'

# Existe uma relacao  entre o valor p de 'Z' e o intervalo de confianca,
# sendo que o valor p sera menor que 0,05 (se e somente se) o intervalo de confianca nao incluir o valor nulo.
# ou seja, se a media e o intervalo de confianca nao passar pelo valor zero.

# calculo da estatistica 'Z' (Borenstein et al. 2009).  
metadata$Z <- with(metadata, d/SE.d) 
  
############################################################################################################################
# AJUSTE DOS PARAMETROS PARA A LEGENDA DO LADO ESQUERDO DO GRAFICO

# agrupa a coluna 'Author' com a coluna 'Year', separado por virgula
metadata$AuthorYear <- apply( metadata[ , 1:2 ] , 1 , paste , collapse = ", " )
  
# calcula 'Ntotal' das amostras de cada estudo
metadata$Ntotal <- apply(metadata[,3:4],1, sum)
  
# coloca os valores da coluna 'Ntotal' entre parenteses
metadata$Ntotal.parenteres <- with(metadata, paste0("(",Ntotal,")"))

# agrupa em uma coluna, a coluna com 'AuthorYear' e 'Ntotal.parenteses'
metadata$left.legend <- apply (metadata[, c(17,19)], 1, paste, collapse = " ")

############################################################################################################################
# AJUSTE DOS PARAMETROS PARA A LEGENDA DO LADO DIREITO DO GRAFICO

# arredonda o valor da media padronizada de 'd' para duas casas decimais
metadata$d.round <- round(metadata$d,2)
 
# arredoda valor do CI.lower
metadata$CI.lower.round <- round(metadata$CI.lower,2)
  
# arredonda valor do CI.upper
metadata$CI.upper.round <- round(metadata$CI.upper,2)

# agrupa a coluna com os valores da media padronizada com as colunas do intervalo de confianca  
metadata$d.CI <- with(metadata, paste0(d.round, " [", CI.lower.round, "; ", CI.upper.round, "] "))
  
############################################################################################################################
# CRIA GRAFICO META ANALISE

# este algoritmo foi inspirado na aula 5b. Graficos avancados
# http://ecologia.ib.usp.br/bie5782/doku.php?id=bie5782:03_apostila:10-graficos02
# porem, foram realizadas varias mudancas para que o grafico ficasse adequado com a proposta da funcao

# cria uma matriz que sera a area do grafico
layout(matrix(c(1,2), ncol=2, nrow=1), width=c(8, 2))

# mostra matriz criada
layout.show(2)

# estabelece os valores da margem do grafico
par (mar=c(5, 10 ,4, 4)) 
  
# plota area do grafico, com os limites de xlim e ylim determinados pelos valores do intervalo de confinca e os valores 
# da media, e ajusta alguns parametros da legenda do grafico
plot(x=NULL, y=NULL, xlim=c(min(metadata$CI.lower -1), max(metadata$CI.upper + 1)), ylim=c(1, nrow(metadata)), type="n", yaxt="n",
     ylab="", xlab="Standardised mean difference (95% CI)",bty= "n", main = "Author, year (total sample)                                                                              SMD [95% CI]", cex=1)

# cria legenda do lado direito (side 2), com os nomes dos autores de cada estudo, ano e o N amostral total (controle + tratamento) para cada estudo
axis(side=2, at=rev(c(1:nrow(metadata))), labels=c(metadata$left.legend), las=2, lwd = 0)

# cria legenda do lado esquerdo (side 4), com os valores da media padronizada, calculada com a diferenca entre os estudos
# e o intervalo de confianca de 95% calculado, com minimo e maximo
axis(side=4, at=rev(c(1:nrow(metadata))), labels=c(metadata$d.CI), las=2, lwd = 0)

# cria linha pontilhada vertical na posicao zero do eixo 'x' 
abline (v=0, lty=2) 

# coloca ponto no valor que foi calculado a media padronizada da diferenca entre os estudos  
points (x=c(metadata$d), y=rev(c(1:nrow(metadata))), pch=19 , cex=2) 

# coloca uma barra vertical no minimo e maximo dos intervalos de confianca  
points (x=c(metadata$CI.lower, metadata$CI.upper), y=rep(rev(c(1:nrow(metadata))), times=2), pch= "|" )
  
# cria segmentos com o minimo e maximo do intervalo de confianca
segments(x0=c(metadata$CI.lower), y0=rev(c(1:nrow(metadata))),x1=c(metadata$CI.upper),y1=rev(c(1:nrow(metadata))))
  
# salva grafico estilo forest plot
png(filename = "metaR_forestplot.png",width = 900, height = 900, units = "px", pointsize = 12,bg = "transparent")

# fim, citacao aula 5b. Graficos avancados

# seleciona apenas as colunas para o data.frame de saida para o usuario
metadata.print <- with(metadata,metadata[,1:16])

# salva na area de trabalho do usuario o forest plot gerado pela funcao
write.table(metadata.print, "metaR_data.csv",sep = ";",col.names = NA, row.names = TRUE)

# salva tabela de saida para usuario com as colunas adicionais que foram realizadas os calculos estatisticos para 
# analise dos estudos na meta-analise

# retorna para o usuario a tabela com os dados que foram inseridos e adicionados pelos calculos pela estatistica de interesse
return(metadata.print)

}

# FIM FUNCAO

###########################################################################################################################
# REFERENCIAS UTILIZADAS PARA O DESENVOLVIMENTO DA FUNCAO

# Adalardo, A. O, e colaboradores. (2018).http://ecologia.ib.usp.br/bie5782/doku.php?id=start
# Crawley, M.J. (2007). The R Book. Imperial College London at Silwood Park, UK
# Curtis, P.S. (1996). A meta-analysis of leaf gas exchange and nitrogen in trees grown under elevated carbon dioxide. Plant, Cell and Environment
# Hedges, L.V., Olkin, I. (1985). Statistical Methods for Meta-Analysis.United Kingdom Edition 
# Koricheva, J., Gurevitch, H., Mengersen. (2013). Handbook of Meta-analysis in ecology and evolution. Princeton University Press.

# Além, das duvidas tiradas no forum, leitura do material disponibilizado e 
# conversas com os colegas da disciplina.
# http://ecologia.ib.usp.br/bie5782/doku.php?id=start