# 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