Meta análise Nutrição animal (Abordagem R do Artigo St Pierre (2001))

meta-analiseEm 2001 N.R St-Pierre publicou uma revisão a convite do periódico Dairy Science intitulada “Integrating Quantitative Findings from Multiple Studies Using Mixed Model Methodology” onde discorreu sobre  o uso de meta-análise para integrar achados científicos. O artigo foi muito bem aceito pela comunidade científica, segundo dados do site scopo preview, houve 322 citações desse artigo. No Brasil, alguns artigos elaborados foram baseados nessas recomendações (Azevedo et al,2010,Vieira et al,2013, Souza,2013).

Embora a meta análise seja algo muito mais complexo que o expresso no artigo e os comandos apresentados não devem ser utilizados roboticamente, ou seja, como uma receita de bolo, esse artigo possui um elevado valor didático para quem quer iniciar seus estudos sobre a técnica de meta anáise.

Acessei esse artigo em 2011, naquela época queria refazer os comandos do artigo no software R, pois estava iniciando o interesse pela meta-análise. No início tive dificuldades, pois não chegava aos mesmos resultados do artigo o que esfriou o projeto de refazê-lo. Logo, descobri que o banco de dados que o autor apresenta ao final do artigo não era o mesmo utilizado para gerar os resultados. Agradeço, ao amigo Luigi Cavalcanti por ter conseguido o banco de dados correto com o pŕopio St Pierre.

Recentemente, encontrei o R-script que havia feito na época e  decidi finalizá-lo e compartilhá-lo aqui no blog, na esperança que aqueles que se interessam também pela meta-análise possam iniciar seus estudos utilizando o R.

Há muitas mais coisas  para se discutir sobre meta análise além do que está no artigo em questão. Assim, pretendo em posts futuros discutí-los no blog.

No post de hoje mostrarei os principais comandos do artigo e as peculiaridades do R para se chegar aos resultados.

Até o próximo post!

#--------------------Pacotes utilizados-------------------------
#  install.packages(c("nlme","car","multcomp","lattice"))
library(nlme)
library(car)
library(multcomp)
library(lattice)
library(lmmfit) #Este arquivo foi removido do repositorio CRAN. Para instalá-lo baixe a versão
#antiga neste endereço: "https://cran.r-project.org/src/contrib/Archive/lmmfit/" e instale-o
#No linux: install.packages("~/Downloads/lmmfit_1.0.tar.gz", repos = NULL, type = "source") ou
# No windows: install.packages("c:\users\nome do usuário\Downloads/lmmfit_1.0.tar.gz", repos = NULL, type = "source")
#Substitua "nome do usuário" pelo nome de usuário do seu sistema.
#----------------------------------------------------------------------------------
DATA<-read.csv("https://www.dropbox.com/s/hibr0cfu91u6xxz/St_Pierre.csv?raw=1")
dados<-groupedData(Y~X|Study,data=DATA)
Lista<-lmList(Y~X,data=dados)
Lista # lista os parametros das regressões
#-------------------------------------------------------------------------------------------#
xyplot(Y~X,DATA) # Gráfico 2a
xyplot(Y~X, groups=Study,type="l",auto.key=list(space="right",title="Estudo"),DATA) #Grafico 2b
#################################################################################
# Modelo que não leva em consideração o efeito dos estudo
#################################################################################
#-------------------------------Figura 3-----------------------------------------
RegSimples<-gls(Y~X,data=dados)#considera os efeitos fixos sem considerar o efeito do estudo
summary(RegSimples)
anova(RegSimples,type="marginal",adjustSigma=F) # Soma de quadrado tipo III
anova(RegSimples,type="sequential",adjustSigma=F) # Soma de quadrado tipo I (sequencial)

#---------------------------------------------------------------------------------------
#Alguns gráficos úteis para avaliar o modelo
#---------------------------------------------------------------------------------------
#Residuos por estudo não estão centrados em zero
plot(RegSimples,Study~resid(.),abline=0,id=0.05,adj=-0.3)
# Gráfico útil para observar onde se encontram os efeitos aleatórios.
# Efeitos aleatórios associados ao intercepto e inclinação (intervalos não se sobrepoem)
plot(intervals(Lista2))

#------------------------------R2 modelo gls ---------------------------------------
#---------------------Função permite estimar o R2 para objetos gls------------------
R2<-function(gnls.obj){
da<-eval(getData(gnls.obj)[,1:2])
resp.name<-all.vars(formula(gnls.obj))
form<-paste(resp.name, "~1", sep="")
m0 <-lm(form, da)
sqn <-sum(resid(gnls.obj)^2)
sqe <-deviance(m0)
r2 <-1-(sqn/sqe)
return(list(R2=r2))
}

R2(RegSimples)
#--------------------------------figura 4 e 5--------------------------------------
plot(RegSimples,Y~X,xlim=c(0,12),ylim=c(-4,20),key=list(corner=c(0.10,0.90),text=list(c("y=a+bx (a=0,b=1)","Y=-5.187620 + 1.959493.X","R²=0.79")),lines=list(lty=c(2,1,0),col=c(1,2,0))),
panel = function(x, y,...) { # Figura 4
panel.xyplot(x,y, ...)
panel.lmline(x,y,col=2,...)
panel.abline(a=c(0.0,1.0),lty=2)
})

plot(RegSimples,resid(.)~X,abline=0,ylab="Predito-Observado",xlab="Predito")#Figura 5
#############################################################################
# Modelo considera o efeito do estudo fixo
################################################################################
#Por padrão o contraste do R é o contr.treatment o qual compara os níveis dos tratamentos
#com o nível base (o primeiro nível). Ele difere do tipo de contraste utilizado
#pelo software SAS, o qual utiliza como base de comparação o último nível ( e não o primeiro
#como ocorre no R). Entender o tipo de contraste utilizado por cada software é importante pois
# impede erros de interpretação dos resultados. Por questões didáticas, modificaremos o tipo
#de contraste padrão do R (função options) para fatores ordenados para "contr.SAS" (o padrão é contr.poly).
#=========================Análise e Resultado (Figura 6)================================
options(contrasts=c(factor="contr.SAS",ordered="contr.SAS"))
RegFixo<-gls(Y~X*Study,data=dados)
(resultado<-summary(RegFixo,verbose=TRUE))
(anovaIII<-anova(RegFixo,type="marginal",adjustSigma=FALSE))
(anovaI<-anova(RegFixo))
R2(RegFixo)
#---------------------------------------------------------------------------------
#-------Obter estimativas do intercepto e inclinação média por estudo--------
#--------------------------------------------------------------------------------
m0 <-gls(Y~0+Study/X,dados)
(LSMEANS<-summary(m0)) # Estimativas de inclinações (efeito de X) para cada estudo
## Matriz de pesos.
#--------------------------------------------------------------
# Estima o intercepto global e testa se seu valor é igual a zero
m<-rbind(rep(1, nlevels(dados$Study))/nlevels(dados$Study))
## Intercepto médio.
X <-cbind(m, 0*m)
GlobIntercept<-summary(glht(m0, linfct=X,rhs=c(0)), test=adjusted(type="none"))
#-------------------------------------------------------------
## Estima a inclinação global e testa se o valor é igual a 1
X1&lt;-cbind(0*m, m)
GlobSlope<-summary(glht(m0, linfct=X1,rhs=c(1)), test=adjusted(type="none")) #Testa se a inclinação do modelo
# é igual a zero.
#---------------------Resultados da figura 6-------------------------------------------------
(Figura6<-list(AnovaI=anovaI,AnovaIII=anovaIII,Resultado=resultado,Médias=LSMEANS,Intercepto_Global=OverIntercpt,Inclinação_Global=OverSlope))
#==================================================================================================
#-------------------------------Não tem no artigo--------------------------------------------------
# Compare este gráfico com o gráfico do modelo anterior. Veja como os resíduos estão mais centrados em zero
plot(RegFixo,Study~resid(.),abline=0,id=0.05,adj=-0.3)
plot(RegFixo,resid(.)~X,abline=0,ylab="Predito-Observado",xlab="Predito")
#--- Compare os modelos pelo AIC (menor melhor), modelo RegFixo preferível-----
anova(RegSimples,RegFixo)
#=====================================fim do resultado (figura 6)===============================
#===============================================================================================
#################################################################################
# MODELO CONSIDERA O EFEITO DO ESTUDO ALEATÓRIO
#################################################################################

# Método gráfico utilizado para determinar qual parâmetro do modelo deve ser incluído no modelo
#aleatório. Repare que os intervalos de confiança obtido para cada estudo não se
#sobrepoem para o intercepto e inclinação, indicando que tanto o intercepto quanto a inclinação devem ser incluídos no modelo aleatório
plot(intervals(Lista))
#---------------------Resultado (Figura 7)--------------------------------
RegMisto<-lme(fixed=Y~1+X,control=controle,random=~X|Study,data=dados,method="REML")
summary(RegMisto,verbose=T)
anova(RegMisto,type="marginal") #anova tipo II
R2(RegMisto)
#Figura 10
plot(RegMisto,fitted(.)~X,id=0.05,adj=0.3,key=list(corner=c(0.10,0.90),text=list(c("Y=-0.61(SE=0.56) + 1.097(SE=0.056).X","R²=0.99"))),
panel = function(x, y,...) {
panel.xyplot(x,y, ...)
panel.lmline(x,y,col=2,...)
})

#--------------------------Avaliação dos resíduos-----------------------------------------
plot(RegMisto,Study~resid(.),abline=0,id=0.05,adj=-0.3)
plot(RegMisto,resid(.)~X,abline=0,ylab="Predito-Observado",xlab="Predito")
#---------------Comparação entre os modelos------------------------------------------------
anova(RegSimples,RegFixo,RegMisto)
#Ops! Ao que parece para esses dados considerar efeito fixo do estudo leva a um
#modelo melhor que o modelo de efeito misto. Este com certeza não era o resultado esperado por St.Pierre
#mas como ele não compara os diferentes modelos em seu artigo, tal fato deve ter passado despercebido no momento
#em que os dados sintéticos foram gerados.Esse fato no entanto não prejudica o valor didático do artigo
#----------------------------------------------------------------------------
# No tópico "Expansion e Reducion" pg 749, o autor fala sobre situações
#onde é possível componentes de variancia e não de covariância serem
#significantes. Nesta situação podemos elaborar um modelo mais simples
#o qual considera a covariância igual a zero. O pacote lme não possui
#um teste de covariância dos parâmetros conforme realizado pelo autor
#utilizando o SAS. Entretanto o pacote lmmfit através da função
#structStep() permite escolher a melhor estrutura da matriz de
#variância-covariância atráves do critério AIC para modelos com um nível
#de agrupamento.

structStep(RegMisto) # pdDiag selecionando, idicando ser os efeitos aleatórios independentes. Esta estrutura é equivalente a estrutura VC utilizada pelo SAS o qual considera diferentes variâncias entre os parâmetros e covariância igual a zero
RegMistoDiag<-lme(fixed=Y~1+X,control=,na.action=na.omit,random=pdDiag(form=~1+X),data=dados,method="REML")
summary(RegMistoDiag,verbose=TRUE)
anova(RegMisto,RegMistoWeight2) # embora o teste L.ration não seja estatisticamente diferentes o modelo 2 é preferível por ter menor AIC

Anúncios
Esse post foi publicado em Dicas R, Meta-Análise, Regressões. Bookmark o link permanente.

Deixe um comentário

Preencha os seus dados abaixo ou clique em um ícone para log in:

Logotipo do WordPress.com

Você está comentando utilizando sua conta WordPress.com. Sair / Alterar )

Imagem do Twitter

Você está comentando utilizando sua conta Twitter. Sair / Alterar )

Foto do Facebook

Você está comentando utilizando sua conta Facebook. Sair / Alterar )

Foto do Google+

Você está comentando utilizando sua conta Google+. Sair / Alterar )

Conectando a %s