[R-br] Desdobrar interação dupla e tripla com glm

Maurício Lordêlo mslordelo em gmail.com
Sexta Outubro 31 20:50:33 BRST 2014


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 em 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 em 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.
>
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20141031/3becfc54/attachment.html>


Mais detalhes sobre a lista de discussão R-br