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

Quinta Novembro 6 22:58:36 BRST 2014

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 em 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 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/20141106/849f6056/attachment.html>