Desdobrar interação dupla e tripla com glm
Saudações a todos! No mês passado, solicitei um auxílio sobre como analisar os dados de um experimento realizado com explantes de bananeiras considerando três fatores. Na ocasião descrevi o problema mais ou menos da seguinte forma (ocorreram algumas mudanças): "*Experimento realizado com explantes de bananeiras considerando três fatores:* *1. Genótipo, com dois níveis;* *2. Meio, com sete níveis;* *3. Tipo de corte, com oito níveis.* *O objetivo é verificar qual tipo de corte induziu melhor o explante na formação do embrião considerando os setes meios e os dois genótipos.* *Para isso foram preparadas duas placas (caracterizando o número de repetições) com uma quantidade de explantes idênticos de acordo com o corte e o genótipo (esta quantidade de explantes em cada placa variou de 6 a 11). Cada placa com estes explantes foi introduzida nos sete diferentes meios. As respostas foram obtidas registrando-se, em relação ao total de cada placa, quantos explantes oxidaram, quantos formaram calos, quantos não desenvolveram, etc.* *Os dados com uma variável resposta está neste link:* *https://www.dropbox.com/s/ci4u8c4vn4cmk4l/dados.txt?dl=0 <https://www.dropbox.com/s/ci4u8c4vn4cmk4l/dados.txt?dl=0> *" Obtive retorno do Walmes que contribui com valiosas sugestões. *O que preciso agora é realizar o desdobramento de interações duplas e triplas utilizando "glm"*. No CRM abaixo, apenas a interação dupla foi significativa. Porém, tenho outras variáveis neste experimento cuja interação tripla foi significativa e será necessário também desdobrar. Na busca que fiz, encontrei apenas comandos de desdobramentos para modelos ajustados com "aov". Antecipadamente, agradeço quem possa contribuir. Maurício dados = read.table("dados.txt",header=TRUE, sep="\t") str(dados) attach(dados) # Modelo fatorial completo, glm (quasi)binomial para variável resposta "n_ex_ox" model1 = glm(cbind(good=n_ex_ox, bad=n_explant-n_ex_ox)~genotipo*meio*tipo, data=dados, family=quasibinomial) anova(model1,test="F") summary(model1)
Minha sugestão ou modelo de trabalhar é considerar funções do pacote doBy com multcomp. Basicamente a doBy::LSmatrix() vai te gerar uma matriz que você vai converter em matriz de contrastes e passar para multcomp::glht(). A glht() admite objetos de classe lm, glm, lme. O detalhe para objetos glm é que as "médias" e comparações são na escala do preditor linear, ou seja, quando você pedir LSmeans() para um binomial não será retornado "média" mas sim o valor esperado de log(1/(1-p)) = X%*%beta, mas aí é só converter em probabilidade com a função inversa da função de ligação. Seguem links que podem servir de inspiração. http://www.leg.ufpr.br/~walmes/analises/MESerafim/abacaxi/abacaxisolo.html http://www.leg.ufpr.br/~walmes/cursoR/geneticaEsalq/script10.html À disposição. Walmes.
Obrigado Walmes! Eu encontrei um post criado por você no "ridículas" e tentei adaptar para o meu problema. Não tenho certeza se está adequado. Acho que não. O motivo é que encontrei interação dupla significativa entre "meio" e "tipo" e no desdobramento final, para identificar os níveis onde há diferença significativa, esta diferença não aparece ( "valor-p" sempre igual a 1). Segue o CMR e agradeço mais uma vez a atenção. *https://www.dropbox.com/s/ci4u8c4vn4cmk4l/dados.txt?dl=0 <https://www.dropbox.com/s/ci4u8c4vn4cmk4l/dados.txt?dl=0>* *require(multcomp)* *require(doBy)* *dados <- read.table("dados.txt",header=TRUE, sep="\t")* *str(dados)* *attach(dados)* *xtabs(~genotipo+meio+tipo, dados)* *## Interação.* *dados$gmt <- interaction(genotipo, meio, tipo)* *attach(dados)* *levels(gmt)* *#Modelo fatorial completo* *m0 = glm(cbind(good=n_ex_ox, bad=n_explant-n_ex_ox)~genotipo*meio*tipo,family="quasibinomial", data=dados)* *par(mfrow=c(2,2))* *plot(m0)* *anova(m0,test="F")* *drop1(m0, test="F", scope=.~.)* *#Modelo só com efeito de primeira grau* *m1 <- update(m0, .~genotipo+meio+tipo)* *anova(m1, m0,test="F")* *## Modelo só com efeitos até segundo grau.* *m2 <- update(m0, .~(genotipo+meio+tipo)^2)* *anova(m2, test="F")* *anova(m2, m0,test="F")* *par(mfrow=c(2,2))* *plot(m2)* *##Obs: Estou estudando o modelo logístico com superdispersão constante para melhorar os resíduos deste modelo. Por enquanto vou utilizar o modelo completo* *do.call(rbind, strsplit(levels(dados$gmt),"\\."))* *M <- LSmatrix(m0, effect=c("tipo","meio"))* *dim(M)* *class(M)* *grid <- attr(M, "grid")* *## Matriz de contrastes * *## Antes executar a Função apc* *apc <- function(pmc, lev=NULL){* * ## apc: all pairwise contrasts, pmc: popMatrix contrast* * nlev <- nrow(pmc)* * rn <- rownames(pmc)* * a <- attr(pmc, "grid")* * if(is.null(lev)){* * if(!is.null(a)){* * lev <- apply(a, 1, paste, collapse=":")* * } else if(!is.null(rn)){* * lev <- rn* * } else {* * lev <- as.character(1:nlev)* * }* * }* * cbn <- combn(seq_along(lev), 2)* * M <- pmc[cbn[1,],]-pmc[cbn[2,],]* * if(is.vector(M)) dim(M) <- c(1, length(M))* * rownames(M) <- paste(lev[cbn[1,]], lev[cbn[2,]], sep="-")* * M* *}* *#Matriz de contrastes entre níveis de `focus` dentro dos níveis de `split`.* *split <- c("meio")* *focus <- "tipo"* *spl <- interaction(grid[,split])* *i <- 1:nrow(grid)* *l <- split(i, f=spl)* *contr <- lapply(l,* * function(row){* * ## Matriz de contrastes par a par.* * a <- apc(M[row,], lev=levels(dados[,focus]))* * ## Prefixo no nome das linhas.* * rownames(a) <- paste(spl[row[1]],* * rownames(a), sep="/")* * return(a)* * })* *contr <- do.call(rbind, contr)* *summary(glht(m0, linfct=contr),test=adjusted(type="fdr"))* Em 30 de outubro de 2014 10:22, walmes . <walmeszeviani@gmail.com> escreveu:
Minha sugestão ou modelo de trabalhar é considerar funções do pacote doBy com multcomp. Basicamente a doBy::LSmatrix() vai te gerar uma matriz que você vai converter em matriz de contrastes e passar para multcomp::glht(). A glht() admite objetos de classe lm, glm, lme. O detalhe para objetos glm é que as "médias" e comparações são na escala do preditor linear, ou seja, quando você pedir LSmeans() para um binomial não será retornado "média" mas sim o valor esperado de log(1/(1-p)) = X%*%beta, mas aí é só converter em probabilidade com a função inversa da função de ligação. Seguem links que podem servir de inspiração.
http://www.leg.ufpr.br/~walmes/analises/MESerafim/abacaxi/abacaxisolo.html http://www.leg.ufpr.br/~walmes/cursoR/geneticaEsalq/script10.html
À disposição. Walmes.
_______________________________________________ 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.
Maurício, Segue código com comentários. ##----------------------------------------------------------------------------- ## Pacotes. require(multcomp) require(doBy) require(lattice) require(plyr) ##----------------------------------------------------------------------------- ## Leitura e reconhecimento dos dados. dados <- read.table("dados.txt", header=TRUE, sep="\t") str(dados) ftable(xtabs(~genotipo+meio+tipo, dados)) dados$prop <- with(dados, n_ex_ox/n_explant) xyplot(prop~tipo|genotipo, groups=meio, data=dados, type=c("p","a")) xyplot(prop~meio|tipo, groups=genotipo, data=dados, type=c("p","a")) ## Elevada presença de celas com proporção amostral igual à zero. Isso ## pode ser problemático pois a proporção amostral (média da amostra) é ## a estimativa de p no caso de uma população sem modelo de ## regressão. Mas p tem que estar no (0,1), não pode ser 0 nem 1. Além ## do mais, o erro padrão da estimativa de p envolve p*(1-p) no ## denominador, então se p -> 0, seu erro padrão tende à Inf. Pode ser ## que a estrutura experimental permita estimar valores não na borda, ## mas isso não é garantido. ##----------------------------------------------------------------------------- ## Ajuste de modelos. ## Modelo fatorial completo (maximal). m0 <- glm(cbind(good=n_ex_ox, bad=n_explant-n_ex_ox)~ genotipo*meio*tipo, family="quasibinomial", data=dados) ## Diagnóstico. par(mfrow=c(2,2)); plot(m0); layout(1) anova(m0,test="F") drop1(m0, test="F", scope=.~.) ## Indicativo de que a tripla pode ser abandonada. ## Modelo só com efeito de primeiro grau. m1 <- update(m0, .~genotipo+meio+tipo) anova(m1, m0,test="F") ## Modelo só com efeitos até segundo grau. m2 <- update(m0, .~(genotipo+meio+tipo)^2) anova(m2, m0, test="F") anova(m2, test="F") ## Indica que a dupla genotipo:tipo pode ser abandonada. m3 <- update(m0, .~meio*(genotipo+tipo)) anova(m3, m0, test="F") anova(m3, m2, test="F") anova(m3, test="F") par(mfrow=c(2,2)); plot(m2); layout(1) summary(m2) summary(m2)$dispersion ##----------------------------------------------------------------------------- ## Comparações múltiplas. ## Usando o modelo m2. M <- LSmatrix(m2, effect=c("tipo", "meio")) ## Pacote wzRfun contém essas funções por conveniência. ## Informação para instalação em https://github.com/walmes/wzRfun. require(wzRfun) grid <- equallevels(attr(M, "grid"), dados) str(grid) ## Fazer comparação de "medias" de tipo estratificado por níveis de ## meio. X <- by(M, INDICES=grid$meio, FUN=as.matrix) X <- lapply(X, "rownames<-", levels(grid$tipo)) str(X) ## lapply(X, apmc, model=m2, focus="tipo", test="none") L <- lapply(X, apmc, model=m2, focus="tipo", test="fdr") L L <- ldply(L) names(L)[1] <- "meio" L <- equallevels(L, dados) str(L) ## names(m2) ## names(m2$family) ##----------------------------------------------------------------------------- ## Problemas. L$p <- m2$family$linkinv(L$estimate) xyplot(p~tipo, groups=meio, data=L, type="b") xyplot(p~meio, groups=tipo, data=L, type="b") plot(sqrt(diag(vcov(m2)))) summary(m2)$coefficients ## Verifica-se muito valor t próximo de zero porque o erro-padrão é ## demasiado grande. ## Por que o teste de comparações múltiplas não rejeitou nenhum ## contraste? Primeiro, o teste F é sobre uma hipótese e as comparações ## múltiplas são para uma série de hipóteses. Não há justificativa para ## que ambos concordem sempre. Nesse caso em particular, parte do ## resultado é porque a medida que estimativa de p vai para 0 ou 1 ## (borda), na escala do preditor corresponde à eta indo para -Inf e Inf ## e o seu erro-padrão se torna cada vez maior. Um contraste tem ## erro-padrão que é função da matriz de covariância das estimativas e ## erros-apdrões altos nela implicam em erros padrões altos para os ## contrastes. Esse tipo de inferência, (assintótica | baseada na matriz ## de covariância | na aproximação quadrática da função objetivo) é ## chamada de Wald. Já a inferência baseada na avaliação da função ## objetivo (que seria a função de log-verossimilhança caso não fosse um ## modelo quasi) não sofre desse problema pois testa hipóteses sem ## depender da matriz de covariância. ## Esse resultado mostra que o planejamento de experimentos com ## respostas do tipo binomial deve antecipar algumas coisas para ser bem ## sucedido. Se espera-se que p esteja próximo das bordas, mais ## repetições devem ser consideradas e também um maior número de ensaios ## Bernoulli. Nesse experimento foram consideradas apenas 2 repetições e ## um número máximo de ensaios de 11. A chance de um p amostral de 0 em ## 11 ensaios é bem maior do que se fossem feitos 50. Talvez, se o ## experimento tivesse 5 repetições com 25 explantes por unidade ## experimental, menor a chance de tais problemas. Sabe-se que elevar o ## número de unidades aumenta o custo e tempo do experimento, porém, se ## ele não atender requisitos mínimos, custo e tempo serão demandados e ## as perguntas do experimentos não serão respondidas. ## Pelos gráficos, tem-se uma sugestão de que o meio M1 difere dos ## demais e que estes tem pouca/nenhuma diferença entre si. Não é ## recomendado criar hipóteses com a análise dos dados e sim antes da ## coleta dos mesmos, mas será explorado fazer fusão de níveis até ## chegar em um modelo reduzído satisfatório já que o procedimento de ## comparações multiplas não foi bem sucedido. ##----------------------------------------------------------------------------- ## Modelos reduzidos pela fusão de níveis dos fatores. dados$meio2 <- dados$meio levels(dados$meio2) levels(dados$meio2) <- c("M5:M7", "M4", "M6", "M2:M3", "M1")[c(5,4,4,2,1,3,1)] xtabs(~meio+meio2, dados) m4 <- glm(cbind(good=n_ex_ox, bad=n_explant-n_ex_ox)~ genotipo*meio2*tipo, family="quasibinomial", data=dados) anova(m0, m4, test="F") anova(m4, test="F") dados$tipo2 <- dados$tipo levels(dados$tipo2) levels(dados$tipo2) <- c("T7", "T5", "T8:T6", "T3", "T2:T4", "T1" )[c(6,5,4,5,2,3,1,3)] xtabs(~tipo+tipo2, dados) m5 <- glm(cbind(good=n_ex_ox, bad=n_explant-n_ex_ox)~ genotipo*meio2*tipo2, family="quasibinomial", data=dados) anova(m0, m5, test="F") anova(m5, test="F") ##----------------------------------------------------------------------------- ## Visualizando as estimativas. M <- LSmatrix(m5, effect=c("tipo2", "meio2", "genotipo")) grid <- equallevels(attr(M, "grid"), dados) str(grid) str(M) ci <- confint(glht(m5, linfct=M), calpha=univariate_calpha()) grid <- cbind(grid, ci$confint) str(grid) grid$p <- m5$family$linkinv(grid$Estimate) xyplot(p~tipo2|genotipo, groups=meio2, data=grid, type="b") xyplot(p~meio2|genotipo, groups=tipo2, data=grid, type="b") xyplot(p~tipo2|meio2, groups=genotipo, data=grid, type="b") À disposição. Walmes.
Muito bom Walmes! Infelizmente o experimento foi desenvolvido com grandes falhas que poderiam ter sido evitadas se os profissionais que atuam na área de Estatística fossem procurados (também) na fase inicial do trabalho. Algumas mudanças deverão ocorrer para o desenvolvimento da análise. Agradeço mais uma vez por ter contribuído com estes valiosos comentários. Maurício Em 1 de novembro de 2014 16:23, walmes . <walmeszeviani@gmail.com> escreveu:
Maurício,
Segue código com comentários.
##----------------------------------------------------------------------------- ## Pacotes.
require(multcomp) require(doBy) require(lattice) require(plyr)
##----------------------------------------------------------------------------- ## Leitura e reconhecimento dos dados.
dados <- read.table("dados.txt", header=TRUE, sep="\t") str(dados)
ftable(xtabs(~genotipo+meio+tipo, dados))
dados$prop <- with(dados, n_ex_ox/n_explant)
xyplot(prop~tipo|genotipo, groups=meio, data=dados, type=c("p","a"))
xyplot(prop~meio|tipo, groups=genotipo, data=dados, type=c("p","a"))
## Elevada presença de celas com proporção amostral igual à zero. Isso ## pode ser problemático pois a proporção amostral (média da amostra) é ## a estimativa de p no caso de uma população sem modelo de ## regressão. Mas p tem que estar no (0,1), não pode ser 0 nem 1. Além ## do mais, o erro padrão da estimativa de p envolve p*(1-p) no ## denominador, então se p -> 0, seu erro padrão tende à Inf. Pode ser ## que a estrutura experimental permita estimar valores não na borda, ## mas isso não é garantido.
##----------------------------------------------------------------------------- ## Ajuste de modelos.
## Modelo fatorial completo (maximal). m0 <- glm(cbind(good=n_ex_ox, bad=n_explant-n_ex_ox)~ genotipo*meio*tipo, family="quasibinomial", data=dados)
## Diagnóstico. par(mfrow=c(2,2)); plot(m0); layout(1)
anova(m0,test="F") drop1(m0, test="F", scope=.~.)
## Indicativo de que a tripla pode ser abandonada.
## Modelo só com efeito de primeiro grau. m1 <- update(m0, .~genotipo+meio+tipo) anova(m1, m0,test="F")
## Modelo só com efeitos até segundo grau. m2 <- update(m0, .~(genotipo+meio+tipo)^2) anova(m2, m0, test="F") anova(m2, test="F")
## Indica que a dupla genotipo:tipo pode ser abandonada. m3 <- update(m0, .~meio*(genotipo+tipo)) anova(m3, m0, test="F") anova(m3, m2, test="F") anova(m3, test="F")
par(mfrow=c(2,2)); plot(m2); layout(1) summary(m2) summary(m2)$dispersion
##----------------------------------------------------------------------------- ## Comparações múltiplas.
## Usando o modelo m2. M <- LSmatrix(m2, effect=c("tipo", "meio"))
## Pacote wzRfun contém essas funções por conveniência. ## Informação para instalação em https://github.com/walmes/wzRfun. require(wzRfun)
grid <- equallevels(attr(M, "grid"), dados) str(grid)
## Fazer comparação de "medias" de tipo estratificado por níveis de ## meio.
X <- by(M, INDICES=grid$meio, FUN=as.matrix) X <- lapply(X, "rownames<-", levels(grid$tipo)) str(X)
## lapply(X, apmc, model=m2, focus="tipo", test="none") L <- lapply(X, apmc, model=m2, focus="tipo", test="fdr") L
L <- ldply(L) names(L)[1] <- "meio" L <- equallevels(L, dados) str(L)
## names(m2) ## names(m2$family)
##----------------------------------------------------------------------------- ## Problemas.
L$p <- m2$family$linkinv(L$estimate) xyplot(p~tipo, groups=meio, data=L, type="b") xyplot(p~meio, groups=tipo, data=L, type="b")
plot(sqrt(diag(vcov(m2)))) summary(m2)$coefficients
## Verifica-se muito valor t próximo de zero porque o erro-padrão é ## demasiado grande.
## Por que o teste de comparações múltiplas não rejeitou nenhum ## contraste? Primeiro, o teste F é sobre uma hipótese e as comparações ## múltiplas são para uma série de hipóteses. Não há justificativa para ## que ambos concordem sempre. Nesse caso em particular, parte do ## resultado é porque a medida que estimativa de p vai para 0 ou 1 ## (borda), na escala do preditor corresponde à eta indo para -Inf e Inf ## e o seu erro-padrão se torna cada vez maior. Um contraste tem ## erro-padrão que é função da matriz de covariância das estimativas e ## erros-apdrões altos nela implicam em erros padrões altos para os ## contrastes. Esse tipo de inferência, (assintótica | baseada na matriz ## de covariância | na aproximação quadrática da função objetivo) é ## chamada de Wald. Já a inferência baseada na avaliação da função ## objetivo (que seria a função de log-verossimilhança caso não fosse um ## modelo quasi) não sofre desse problema pois testa hipóteses sem ## depender da matriz de covariância.
## Esse resultado mostra que o planejamento de experimentos com ## respostas do tipo binomial deve antecipar algumas coisas para ser bem ## sucedido. Se espera-se que p esteja próximo das bordas, mais ## repetições devem ser consideradas e também um maior número de ensaios ## Bernoulli. Nesse experimento foram consideradas apenas 2 repetições e ## um número máximo de ensaios de 11. A chance de um p amostral de 0 em ## 11 ensaios é bem maior do que se fossem feitos 50. Talvez, se o ## experimento tivesse 5 repetições com 25 explantes por unidade ## experimental, menor a chance de tais problemas. Sabe-se que elevar o ## número de unidades aumenta o custo e tempo do experimento, porém, se ## ele não atender requisitos mínimos, custo e tempo serão demandados e ## as perguntas do experimentos não serão respondidas.
## Pelos gráficos, tem-se uma sugestão de que o meio M1 difere dos ## demais e que estes tem pouca/nenhuma diferença entre si. Não é ## recomendado criar hipóteses com a análise dos dados e sim antes da ## coleta dos mesmos, mas será explorado fazer fusão de níveis até ## chegar em um modelo reduzído satisfatório já que o procedimento de ## comparações multiplas não foi bem sucedido.
##----------------------------------------------------------------------------- ## Modelos reduzidos pela fusão de níveis dos fatores.
dados$meio2 <- dados$meio levels(dados$meio2) levels(dados$meio2) <- c("M5:M7", "M4", "M6", "M2:M3", "M1")[c(5,4,4,2,1,3,1)] xtabs(~meio+meio2, dados)
m4 <- glm(cbind(good=n_ex_ox, bad=n_explant-n_ex_ox)~ genotipo*meio2*tipo, family="quasibinomial", data=dados) anova(m0, m4, test="F") anova(m4, test="F")
dados$tipo2 <- dados$tipo levels(dados$tipo2) levels(dados$tipo2) <- c("T7", "T5", "T8:T6", "T3", "T2:T4", "T1" )[c(6,5,4,5,2,3,1,3)] xtabs(~tipo+tipo2, dados)
m5 <- glm(cbind(good=n_ex_ox, bad=n_explant-n_ex_ox)~ genotipo*meio2*tipo2, family="quasibinomial", data=dados) anova(m0, m5, test="F") anova(m5, test="F")
##----------------------------------------------------------------------------- ## Visualizando as estimativas.
M <- LSmatrix(m5, effect=c("tipo2", "meio2", "genotipo")) grid <- equallevels(attr(M, "grid"), dados) str(grid)
str(M) ci <- confint(glht(m5, linfct=M), calpha=univariate_calpha()) grid <- cbind(grid, ci$confint) str(grid)
grid$p <- m5$family$linkinv(grid$Estimate) xyplot(p~tipo2|genotipo, groups=meio2, data=grid, type="b") xyplot(p~meio2|genotipo, groups=tipo2, data=grid, type="b") xyplot(p~tipo2|meio2, groups=genotipo, data=grid, type="b")
À disposição. Walmes.
_______________________________________________ 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.
participantes (2)
-
Maurício Lordêlo -
walmes .