[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