Teste de médias em GLMM com covariável!
Boa tarde senhores! Estou a alguns dias pesquisando em como fazer um teste de comparações em modelo com covariáveis e cujos slopes são diferentes significativamente, já que, as comparações devem ser feitas com as médias ajustadas para as covariáveis. Já procurei no nablle e pela internet, mais nada esclarecedor. Se alguém já trabalhou com este tipo de problema e que queira me dar uma dica, principalmente no que diz respeito a algum material que explique filosoficamente e matricialmente como funciona as comparações múltiplas neste casos de covariável e ficarei extremamente grato. Segue um CRM do que eu vos falo. set.seed(1) da <- data.frame(Trat=gl(5,10),covi=sample(rpois(100,30),50,replace=TRUE),vari=sample(rpois(100,60),50,replace=TRUE),animal=seq(1:50)) da library(lme4) ancova <- glmer(vari ~ Trat*covi + (1|animal%in%Trat),data=da,family=poisson) Anova(ancova)#A interação mostra que as diferenças de slopes entre Trat são significativas ancova2 <- update(ancova,~.-Trat:covi) anova(ancova2,ancova) library(multcomp) summary(glht(ancova,linfct=mcp(Trat='Tukey'),interaction_average=TRUE,covariate_average=TRUE))#Acredito que isto é feito sem correção para as covariáveis. Em algumas buscas, encontrei que a função effect da library ``effects'' faz tal ajutes, mais não estou seguro quanto ao seu uso. Por isso, com o entendimento de como funciona isto matricialmente e acho que consigo resolver via 'glht'. \begin{signature} <<>>= Prof. Dr. Ivan Bezerra Allaman Universidade Estadual de Santa Cruz Departamento de Ciências Exatas e Tecnológicas Ilhéus/BA - Brasil Fone: +55 73 3680-5596 E-mail: ivanalaman@yahoo.com.br/ivanalaman@gmail.com @ \end{signature}
Ivan, Vamos primeiro primeiro apodar as arestas a fim de fazer a coisa rolar. Você tem um experimento com um fator de 5 níveis nominais, denominado de F. Cada animal é uma unidade experimental do qual se observa uma resposta y e uma covariável x. Acredita-se que x afeta y e portanto quer-se remover esse efeito para poder estudar a variação explicada por F. Pode-se supor que F e x tenham efeitos aditivos (sem interação) ou pode-se assumir que existe interação. Na ausência de interação as diferenças entre níveis de F são as mesmas sem importar o nível de x, mas o mesmo não vale para presença de interação, de forma que, você deve definir antecipadamente em que nível(is) de x quer comparar níveis de F. Não sei o CMR é exatamente o que você tem, pois o CMR não parece ser um caso de modelo misto porque animal é a unidade experimental e não um fator de agrupamento (embora para Poisson este seja um modelo estimável), enfim... vou considerar um glm fixo, daí é só fazer as devidas extensões para o glmm. da <- data.frame(F=gl(5,20), x=runif(5*20)) da$Ey <- with(da, exp(1+as.numeric(F)/2+(as.numeric(F)-2)/3*x)) da$y <- with(da, rpois(F, lambda=exp(1+as.numeric(F)/2+(as.numeric(F)-2)/3*x))) with(da, tapply(y, F, mean)) require(lattice) xyplot(y~x, groups=F, data=da, type="a") ac0 <- glm(y~F*x, data=da, family=poisson) anova(ac0, test="Chisq") betas <- coef(ac0) require(doBy) # médias populacionais marginais em x=0.5 mpm <- popMatrix(ac0, effect="F", at=list(x=0.5)) mpm%*%betas exp(mpm%*%betas) # diferença nas médias populacionais marginais duas à duas nlevels(da$F) comp <- t(combn(1:nlevels(da$F), 2)) rownames(comp) <- apply(comp, 1, paste, collapse="-") dmpm <- t(apply(comp, 1, function(vs) mpm[vs[1],]-mpm[vs[2],])) rownames(dmpm) <- rownames(comp) apply(dmpm, 1, function(x) x%*%betas) require(multcomp) summary(glht(ac0, linfct=mpm)) # médias populacionais marginais summary(glht(ac0, linfct=dmpm)) # contraste entre as médias populacionais marginais # na verdade isso não é a média, mas sim o exp() disso E é isso aí Ivan. Esse procedimento aí que fiz foi o que deixou um determinado software famoso por dar esses resultados aos usuários. Trata-se de algo que, depois de fazer pela primeira vez e entender, passa a ser trivial. À disposição. Walmes. ========================================================================== Walmes Marques Zeviani LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W) Departamento de Estatística - Universidade Federal do Paraná fone: (+55) 41 3361 3573 VoIP: (3361 3600) 1053 1173 e-mail: walmes@ufpr.br twitter: @walmeszeviani homepage: http://www.leg.ufpr.br/~walmes linux user number: 531218 ==========================================================================
Walmes, o mestre da didática. Parabéns!! Walmes não deu para assistir a sua comunicação. Esepro que tenha transcorrido tudo bem. Abs Em 9 de maio de 2012 16:05, Walmes Zeviani <walmeszeviani@gmail.com>escreveu:
Ivan,
Vamos primeiro primeiro apodar as arestas a fim de fazer a coisa rolar. Você tem um experimento com um fator de 5 níveis nominais, denominado de F. Cada animal é uma unidade experimental do qual se observa uma resposta y e uma covariável x. Acredita-se que x afeta y e portanto quer-se remover esse efeito para poder estudar a variação explicada por F. Pode-se supor que F e x tenham efeitos aditivos (sem interação) ou pode-se assumir que existe interação. Na ausência de interação as diferenças entre níveis de F são as mesmas sem importar o nível de x, mas o mesmo não vale para presença de interação, de forma que, você deve definir antecipadamente em que nível(is) de x quer comparar níveis de F. Não sei o CMR é exatamente o que você tem, pois o CMR não parece ser um caso de modelo misto porque animal é a unidade experimental e não um fator de agrupamento (embora para Poisson este seja um modelo estimável), enfim... vou considerar um glm fixo, daí é só fazer as devidas extensões para o glmm.
da <- data.frame(F=gl(5,20), x=runif(5*20)) da$Ey <- with(da, exp(1+as.numeric(F)/2+(as.numeric(F)-2)/3*x)) da$y <- with(da, rpois(F, lambda=exp(1+as.numeric(F)/2+(as.numeric(F)-2)/3*x))) with(da, tapply(y, F, mean))
require(lattice) xyplot(y~x, groups=F, data=da, type="a")
ac0 <- glm(y~F*x, data=da, family=poisson) anova(ac0, test="Chisq") betas <- coef(ac0)
require(doBy) # médias populacionais marginais em x=0.5 mpm <- popMatrix(ac0, effect="F", at=list(x=0.5)) mpm%*%betas exp(mpm%*%betas)
# diferença nas médias populacionais marginais duas à duas nlevels(da$F) comp <- t(combn(1:nlevels(da$F), 2)) rownames(comp) <- apply(comp, 1, paste, collapse="-")
dmpm <- t(apply(comp, 1, function(vs) mpm[vs[1],]-mpm[vs[2],])) rownames(dmpm) <- rownames(comp)
apply(dmpm, 1, function(x) x%*%betas)
require(multcomp) summary(glht(ac0, linfct=mpm)) # médias populacionais marginais summary(glht(ac0, linfct=dmpm)) # contraste entre as médias populacionais marginais
# na verdade isso não é a média, mas sim o exp() disso
E é isso aí Ivan. Esse procedimento aí que fiz foi o que deixou um determinado software famoso por dar esses resultados aos usuários. Trata-se de algo que, depois de fazer pela primeira vez e entender, passa a ser trivial.
À disposição. Walmes.
========================================================================== Walmes Marques Zeviani LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W) Departamento de Estatística - Universidade Federal do Paraná fone: (+55) 41 3361 3573 VoIP: (3361 3600) 1053 1173 e-mail: walmes@ufpr.br twitter: @walmeszeviani homepage: http://www.leg.ufpr.br/~walmes linux user number: 531218 ==========================================================================
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível.
-- Guilherme Moraes Ferraudo Jaboticabal/Campinas - SP http://lattes.cnpq.br/2096118558794430
Ferraudo, tudo bem? Cara, quando li a mensagem, pensei a mesma coisa! Cada dia aprendo mais com o Walmes. Parabéns, Walmes. Poucas vezes li uma explicação tão simples e didática de uma análise de covariâncias. Abs Olympio Em 09/05/2012 18:29, Guilherme Moraes Ferraudo escreveu:
Walmes, o mestre da didática. Parabéns!!
Walmes não deu para assistir a sua comunicação. Esepro que tenha transcorrido tudo bem.
Abs
Em 9 de maio de 2012 16:05, Walmes Zeviani <walmeszeviani@gmail.com <mailto:walmeszeviani@gmail.com>> escreveu:
Ivan,
Vamos primeiro primeiro apodar as arestas a fim de fazer a coisa rolar. Você tem um experimento com um fator de 5 níveis nominais, denominado de F. Cada animal é uma unidade experimental do qual se observa uma resposta y e uma covariável x. Acredita-se que x afeta y e portanto quer-se remover esse efeito para poder estudar a variação explicada por F. Pode-se supor que F e x tenham efeitos aditivos (sem interação) ou pode-se assumir que existe interação. Na ausência de interação as diferenças entre níveis de F são as mesmas sem importar o nível de x, mas o mesmo não vale para presença de interação, de forma que, você deve definir antecipadamente em que nível(is) de x quer comparar níveis de F. Não sei o CMR é exatamente o que você tem, pois o CMR não parece ser um caso de modelo misto porque animal é a unidade experimental e não um fator de agrupamento (embora para Poisson este seja um modelo estimável), enfim... vou considerar um glm fixo, daí é só fazer as devidas extensões para o glmm.
da <- data.frame(F=gl(5,20), x=runif(5*20)) da$Ey <- with(da, exp(1+as.numeric(F)/2+(as.numeric(F)-2)/3*x)) da$y <- with(da, rpois(F, lambda=exp(1+as.numeric(F)/2+(as.numeric(F)-2)/3*x))) with(da, tapply(y, F, mean))
require(lattice) xyplot(y~x, groups=F, data=da, type="a")
ac0 <- glm(y~F*x, data=da, family=poisson) anova(ac0, test="Chisq") betas <- coef(ac0)
require(doBy) # médias populacionais marginais em x=0.5 mpm <- popMatrix(ac0, effect="F", at=list(x=0.5)) mpm%*%betas exp(mpm%*%betas)
# diferença nas médias populacionais marginais duas à duas nlevels(da$F) comp <- t(combn(1:nlevels(da$F), 2)) rownames(comp) <- apply(comp, 1, paste, collapse="-")
dmpm <- t(apply(comp, 1, function(vs) mpm[vs[1],]-mpm[vs[2],])) rownames(dmpm) <- rownames(comp)
apply(dmpm, 1, function(x) x%*%betas)
require(multcomp) summary(glht(ac0, linfct=mpm)) # médias populacionais marginais summary(glht(ac0, linfct=dmpm)) # contraste entre as médias populacionais marginais
# na verdade isso não é a média, mas sim o exp() disso
E é isso aí Ivan. Esse procedimento aí que fiz foi o que deixou um determinado software famoso por dar esses resultados aos usuários. Trata-se de algo que, depois de fazer pela primeira vez e entender, passa a ser trivial.
À disposição. Walmes.
========================================================================== Walmes Marques Zeviani LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W) Departamento de Estatística - Universidade Federal do Paraná fone: (+55) 41 3361 3573 <tel:%28%2B55%29%2041%203361%203573> VoIP: (3361 3600) 1053 1173 e-mail: walmes@ufpr.br <mailto:walmes@ufpr.br> twitter: @walmeszeviani homepage: http://www.leg.ufpr.br/~walmes <http://www.leg.ufpr.br/%7Ewalmes> linux user number: 531218 ==========================================================================
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br <mailto:R-br@listas.c3sl.ufpr.br> https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível.
-- Guilherme Moraes Ferraudo Jaboticabal/Campinas - SP http://lattes.cnpq.br/2096118558794430
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível.
Td otimo Olympio. Abs Em 9 de maio de 2012 20:54, Olympio Neto <olympiotneto@gmail.com> escreveu:
Ferraudo, tudo bem?
Cara, quando li a mensagem, pensei a mesma coisa! Cada dia aprendo mais com o Walmes. Parabéns, Walmes. Poucas vezes li uma explicação tão simples e didática de uma análise de covariâncias.
Abs
Olympio
Em 09/05/2012 18:29, Guilherme Moraes Ferraudo escreveu:
Walmes, o mestre da didática. Parabéns!!
Walmes não deu para assistir a sua comunicação. Esepro que tenha transcorrido tudo bem.
Abs
Em 9 de maio de 2012 16:05, Walmes Zeviani <walmeszeviani@gmail.com>escreveu:
Ivan,
Vamos primeiro primeiro apodar as arestas a fim de fazer a coisa rolar. Você tem um experimento com um fator de 5 níveis nominais, denominado de F. Cada animal é uma unidade experimental do qual se observa uma resposta y e uma covariável x. Acredita-se que x afeta y e portanto quer-se remover esse efeito para poder estudar a variação explicada por F. Pode-se supor que F e x tenham efeitos aditivos (sem interação) ou pode-se assumir que existe interação. Na ausência de interação as diferenças entre níveis de F são as mesmas sem importar o nível de x, mas o mesmo não vale para presença de interação, de forma que, você deve definir antecipadamente em que nível(is) de x quer comparar níveis de F. Não sei o CMR é exatamente o que você tem, pois o CMR não parece ser um caso de modelo misto porque animal é a unidade experimental e não um fator de agrupamento (embora para Poisson este seja um modelo estimável), enfim... vou considerar um glm fixo, daí é só fazer as devidas extensões para o glmm.
da <- data.frame(F=gl(5,20), x=runif(5*20)) da$Ey <- with(da, exp(1+as.numeric(F)/2+(as.numeric(F)-2)/3*x)) da$y <- with(da, rpois(F, lambda=exp(1+as.numeric(F)/2+(as.numeric(F)-2)/3*x))) with(da, tapply(y, F, mean))
require(lattice) xyplot(y~x, groups=F, data=da, type="a")
ac0 <- glm(y~F*x, data=da, family=poisson) anova(ac0, test="Chisq") betas <- coef(ac0)
require(doBy) # médias populacionais marginais em x=0.5 mpm <- popMatrix(ac0, effect="F", at=list(x=0.5)) mpm%*%betas exp(mpm%*%betas)
# diferença nas médias populacionais marginais duas à duas nlevels(da$F) comp <- t(combn(1:nlevels(da$F), 2)) rownames(comp) <- apply(comp, 1, paste, collapse="-")
dmpm <- t(apply(comp, 1, function(vs) mpm[vs[1],]-mpm[vs[2],])) rownames(dmpm) <- rownames(comp)
apply(dmpm, 1, function(x) x%*%betas)
require(multcomp) summary(glht(ac0, linfct=mpm)) # médias populacionais marginais summary(glht(ac0, linfct=dmpm)) # contraste entre as médias populacionais marginais
# na verdade isso não é a média, mas sim o exp() disso
E é isso aí Ivan. Esse procedimento aí que fiz foi o que deixou um determinado software famoso por dar esses resultados aos usuários. Trata-se de algo que, depois de fazer pela primeira vez e entender, passa a ser trivial.
À disposição. Walmes.
========================================================================== Walmes Marques Zeviani LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W) Departamento de Estatística - Universidade Federal do Paraná fone: (+55) 41 3361 3573 VoIP: (3361 3600) 1053 1173 e-mail: walmes@ufpr.br twitter: @walmeszeviani homepage: http://www.leg.ufpr.br/~walmes linux user number: 531218 ==========================================================================
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível.
-- Guilherme Moraes Ferraudo Jaboticabal/Campinas - SP http://lattes.cnpq.br/2096118558794430
_______________________________________________ R-br mailing listR-br@listas.c3sl.ufpr.brhttps://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível.
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível.
-- Guilherme Moraes Ferraudo Jaboticabal/Campinas - SP http://lattes.cnpq.br/2096118558794430
Boa noite Walmes, obrigado pela pronta resposta! É provável que não estejamos falando do mesmo tipo de ANCOVA (Análise de covariância). No entanto, como você está anos luz a frente, solucione algumas dúvidas primeiro para que eu não fale besteira. Eu não entendi a seguinte pedaço de frase:``você deve definir antecipadamente em que nível(is) de x quer comparar níveis de F.'' O tipo de ANCOVA que estou falando, não tem sentido comparar níveis de um grupo para um determinado nível de covariável, afinal de contas ela é quantitativa e sua inclusão foi muito bem explicada por você nos trechos acima. O fato de incluirmos a interação entre a covariável e o fator é devido a intenção de avaliar se será necessário ajustar as médias considerandos slopes diferentes (caso dê significativo), ou se considera um slope comum na correção das médias. Até cheguei cogitar de no ítem `` médias populacionais marginais em x=0.5'', você ter considerado x=0.5 como um slope comum para ajuste das médias, mais pelos resultados, creio que não foi essa a intenção. Vamos pegar o gancho do seu CRM, e discutindo ao longo das análises para que eu possa me fazer entender. Considerando já todo o enredo brilhantemente explicado por ti de acordo com o CRM temos o seguinte modelo: ac0 <- glm(y~F*x, data=da, family=poisson) anova(ac0, test="Chisq") De acordo com os resultados, há a necessidade de se ajustar as médias considerando uma regressão diferente para cada grupo. E de fato, os slopes são diferentes estimados para cada grupo. Graficamente, de fato: require(lattice) xyplot(y~x, groups=F, data=da, type="a") Analiticamente, de fato: slopes <- list() for(i in 1:5){ slopes[[i]] <- coef(glm(y ~ x, data=subset(da,F==i)))[2] } slopes Segundo a literatura, o ideal é que ao se comparar as médias, estas devem estar ajustadas para as covariáveis (caso ela for significativa é claro). Então ao invés de eu comparar estas médias: with(da, tapply(y, F, mean)) e devo comparar estas médias: novamedia_F1 <- with(da,mean(y[F==1])) - (exp(as.numeric(slopes[1]))*(with(da,mean(x[F==1]))-mean(da$x))) novamedia_F2 <- with(da,mean(y[F==2])) - (exp(as.numeric(slopes[2]))*(with(da,mean(x[F==2]))-mean(da$x))) novamedia_F3 <- with(da,mean(y[F==3])) - (exp(as.numeric(slopes[3]))*(with(da,mean(x[F==3]))-mean(da$x))) novamedia_F4 <- with(da,mean(y[F==4])) - (exp(as.numeric(slopes[4]))*(with(da,mean(x[F==4]))-mean(da$x))) novamedia_F5 <- with(da,mean(y[F==5])) - (exp(as.numeric(slopes[5]))*(with(da,mean(x[F==5]))-mean(da$x))) rbind(novamedia_F1,novamedia_F2,novamedia_F3,novamedia_F4,novamedia_F5) E a grande pergunta: Como eu posso fazer isso no R? \begin{signature} <<>>= Prof. Dr. Ivan Bezerra Allaman Universidade Estadual de Santa Cruz Departamento de Ciências Exatas e Tecnológicas Ilhéus/BA - Brasil Fone: +55 73 3680-5596 E-mail: ivanalaman@yahoo.com.br/ivanalaman@gmail.com @ \end{signature} ________________________________ De: Walmes Zeviani <walmeszeviani@gmail.com> Para: r-br@listas.c3sl.ufpr.br; Ivan Bezerra Allaman <ivanalaman@yahoo.com.br> Enviadas: Quarta-feira, 9 de Maio de 2012 16:05 Assunto: Re: [R-br] Teste de médias em GLMM com covariável! Ivan, Vamos primeiro primeiro apodar as arestas a fim de fazer a coisa rolar. Você tem um experimento com um fator de 5 níveis nominais, denominado de F. Cada animal é uma unidade experimental do qual se observa uma resposta y e uma covariável x. Acredita-se que x afeta y e portanto quer-se remover esse efeito para poder estudar a variação explicada por F. Pode-se supor que F e x tenham efeitos aditivos (sem interação) ou pode-se assumir que existe interação. Na ausência de interação as diferenças entre níveis de F são as mesmas sem importar o nível de x, mas o mesmo não vale para presença de interação, de forma que, você deve definir antecipadamente em que nível(is) de x quer comparar níveis de F. Não sei o CMR é exatamente o que você tem, pois o CMR não parece ser um caso de modelo misto porque animal é a unidade experimental e não um fator de agrupamento (embora para Poisson este seja um modelo estimável), enfim... vou considerar um glm fixo, daí é só fazer as devidas extensões para o glmm. da <- data.frame(F=gl(5,20), x=runif(5*20)) da$Ey <- with(da, exp(1+as.numeric(F)/2+(as.numeric(F)-2)/3*x)) da$y <- with(da, rpois(F, lambda=exp(1+as.numeric(F)/2+(as.numeric(F)-2)/3*x))) with(da, tapply(y, F, mean)) require(lattice) xyplot(y~x, groups=F, data=da, type="a") ac0 <- glm(y~F*x, data=da, family=poisson) anova(ac0, test="Chisq") betas <- coef(ac0) require(doBy) # médias populacionais marginais em x=0.5 mpm <- popMatrix(ac0, effect="F", at=list(x=0.5)) mpm%*%betas exp(mpm%*%betas) # diferença nas médias populacionais marginais duas à duas nlevels(da$F) comp <- t(combn(1:nlevels(da$F), 2)) rownames(comp) <- apply(comp, 1, paste, collapse="-") dmpm <- t(apply(comp, 1, function(vs) mpm[vs[1],]-mpm[vs[2],])) rownames(dmpm) <- rownames(comp) apply(dmpm, 1, function(x) x%*%betas) require(multcomp) summary(glht(ac0, linfct=mpm)) # médias populacionais marginais summary(glht(ac0, linfct=dmpm)) # contraste entre as médias populacionais marginais # na verdade isso não é a média, mas sim o exp() disso E é isso aí Ivan. Esse procedimento aí que fiz foi o que deixou um determinado software famoso por dar esses resultados aos usuários. Trata-se de algo que, depois de fazer pela primeira vez e entender, passa a ser trivial. À disposição. Walmes. ========================================================================== Walmes Marques Zeviani LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W) Departamento de Estatística - Universidade Federal do Paraná fone: (+55) 41 3361 3573 VoIP: (3361 3600) 1053 1173 e-mail: walmes@ufpr.br twitter: @walmeszeviani homepage: http://www.leg.ufpr.br/~walmes linux user number: 531218 ==========================================================================
participantes (4)
-
Guilherme Moraes Ferraudo -
Ivan Bezerra Allaman -
Olympio Neto -
Walmes Zeviani